Regular Article
Dynamical wellkilling simulation of a vertical H_{2}Scontaining natural gas well
^{1}
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Sichuan, Chengdu 610500, China
^{2}
Petroleum Engineering School, Southwest Petroleum University, Sichuan, Chengdu 610500, China
^{*} Correspondence authors: maoliangjie@qq.com, liuqy66@aliyun.com
Received:
3
March
2020
Accepted:
27
August
2020
This work aims to explore the dynamical wellkilling process of a vertical H_{2}Scontaining natural gas well. A dynamical wellkilling model considering an H_{2}S solubility was established to simulate the overflow and wellkilling process of a vertical H_{2}Scontaining natural gas well. The mass and momentum equations of the coupled model were solved using finite difference method, while the transient temperature prediction model was solved using finite volume method. The coupled model was validated by reproducing experimental data and field data of Well Tiandong #5. The effect of H_{2}S content, mud displacement, drilling fluid density, and initial overflow volume on the dynamical wellkilling process of an H_{2}Scontaining natural gas well were obtained and analyzed in this work. Results showed that H_{2}S will gasify near wellhead during well killing when casing pressure decreases. To balance the bottom hole pressure, when H_{2}S releases, the casing pressure increases as H_{2}S content increases. As initial overflow volume increases, the annular temperature, annular pressure and the casing pressure increase significantly. When H_{2}S gasifies, the casing pressure applied at wellhead should be higher at lower initial overflow volume to balance bottom hole pressure. In the wellkilling process, the annular pressure and temperature decrease as drilling fluid density increases and a lower casing pressure is needed for balancing bottom hole pressure. The casing pressure is lower at a higher displacement for higher friction resistance. Besides, as wellkilling displacement increases H_{2}S will gasify at an earlier time. When drilling for H_{2}Scontaining natural gas well, early detection of gas kick should be more frequent to avoid severe overflow. Besides, higher displacement and density of drilling fluid should be considered to avoid stratum fracturing and prevent leakage accidents under the premise of meeting drilling requirements.
© L. Mao et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
ρ_{l} : Density of drilling fluid
v_{l} : Velocity of drilling fluid
g : Local acceleration of gravity
c_{1} : Specific heat capacity of the drilling fluid in the drilling string
T_{1} : Drilling fluid temperature in the drilling string
u_{p} : Velocity in the x direction of the drilling fluid in the drilling string
v_{p} : Velocity in the y direction of the drilling fluid in the drilling string
Γ_{x}: Heat transfer coefficient of the drilling fluid in the x direction
Γ_{y}: Heat transfer coefficient of the drilling fluid in the y direction
ρ_{2} : Density of the drilling string
T_{2}: Drilling string temperature
V_{pg} : Pit gain at shutin time
v_{sg} : Superficial velocities of gas
v_{sl} : Superficial velocities of liquid
k_{1}, k_{2}: Correction coefficients
y_{i} : Mole fractions of H_{2}S in the gas
x_{i} : Mole fractions of H_{2}S in the liquid
: Fugacity coefficients of H_{2}S in the gas
: Fugacity coefficients of H_{2}S in the liquid
k_{ij} : Interaction coefficient between H_{2}S and hydrocarbon components in natural gas
T_{ci} : Critical temperature of H_{2}S
P_{ci} : Critical pressure of H_{2}S
T_{ri} : Correspondent temperature of H_{2}S
C_{0} : Distribution coefficient
v_{m} : Mixed velocity of gas and liquid
n : Flow index of the drilling fluid
R_{ e } : Reynolds number of the mixed fluids
ε_{e} : Equivalent absolute roughness
W_{g} : Molar mass of natural gas
P_{la} : Circulating pressure loss
P_{ma} : Hydrostatic pressure created by the fluid column
R_{ e } : Mean Reynolds number of the mixed fluids
: Fugacity coefficients of H_{2}S in the gas
: Fugacity coefficients of H_{2}S in the liquid
Z : Natural gas compression factor
1 Introduction
Gas kick and wellkilling are the main concerns of well control during the drilling of oil and gas fields. Various multiphase flow models have been proposed to predict the flow pattern distribution in the wellbore after gas kick (AmayaGómez et al., 2019; Ansari et al., 1994; Bilicki and Kestin, 1987; Ebrahimi and Khamehchi, 2015; Feng et al., 2015; Guo et al., 2018; Hasan and Kabir, 1988; Kelessidis and Dukler, 1989; Li and Mouline, 1997; Zhang et al., 2019); among them, the most widely used in drilling engineering are the patterns of bubble, slug, churn, and annular flows (Osman, 2004; Raghavan, 1989; Sun et al., 2017). On the basis of these studies, many researchers have focused on the multiphase flow in the wellbore under different conditions. Nickens (1987) proposed a dynamic computer model by a solution of the mass and momentum balance equations for predicting the dynamical twophase flow behavior and pressure response of the wellbore. Sun and Wang proposed a mathematical model to predict the multiphase flow in a wellbore for deep water wells (Wang et al., 2016). Their model has also been used to study the phase transition of hydrate during the exploitation of natural gas hydrate (Wang and Sun, 2009, 2014). Xu et al. (2018, 2019) developed a nonisothermal twophase flow model to investigate the effect of major parameters on the twophase flow behavior in the wellbore. They found that temperature, pressure, and solubility fields are mutually influential. The gas solubility effect and heat transfer effect influence gas kick characteristics significantly. Meng et al. (2015) proposed a transient gas–liquid drift flux model based on Shi slip relation for simulating gas kick/injection and solved the model numerically using finite volume method and advection upstream splitting method V scheme. Xie et al. (2014) established a transient model to study the changing rule of wellhead parameters during blowout and found that the gas outlet velocity increases as well depth increases. Hou et al. (2019) conducted an experiment to study the heat transfer for gasliquid twophase flow in wellbore and developed a new model of convective heat transfer coefficient. Fu et al. (2019) proposed the mathematic model of wellbore annulus transient water hammer with the consideration of transient multiphase flow characteristics and pointed out that both the additional water hammer pressure and shutin casing pressure generated by the closure of BOP are likely to cause formation at the shallow casing shoe damage.
With the increasing demand for natural gas, an increasing number of gas reservoirs that contain H_{2}S will be drilled (Guo and Wang, 2016; Yin et al., 2015; Zhang et al., 2020; Zhu et al., 2011). The existing states of H_{2}S are greatly affected by temperature and pressure. The state may transform from supercritical to gas in the wellbore as annular pressure and temperature decrease, and this condition influences annular pressure and flow pattern distribution. Therefore, the typical gas–liquid twophase flow models may be unsuitable for predicting the multiphase flow behavior in H_{2}Scontaining natural gas wells. The “12·23” Gas Blowout Accident in Kaixian County, Chongqing has elicited scholarly attention on phase transition of H_{2}S in well control. Some scholars have focused on multiphase flow behavior of H_{2}Scontaining natural gas wells. Sun and Wang et al. established a multiphase flow model that considers the phase change of the H_{2}S in the drilling fluid; simulation results showed that, as the invasion gas moves up along the wellbore to the critical position of H_{2}S, the released H_{2}S may cause a rapid volume expansion, thereby increasing the blowout risk (Sun et al., 2013, 2018).
It can be seen that most of the existing studies are focused on the gas kick simulation of a well and, the wellkilling simulation after gas kick is scarce, especially for an H_{2}Scontaining natural gas well. By applying casing pressure in the wellhead during well killing, the twophase flow behavior is different from that in the overflow (gas kick). Thus, the phase change behavior of H_{2}S may be different from that described in the abovementioned published literature. Thus, it is of significance to study the wellkilling process of an H_{2}Scontaining natural gas well to offer a guidance for exploitation of H_{2}Scontaining gas reservoir.
Thus, the aim of this study is to establish a dynamical wellkilling model considering an H_{2}S solubility to simulate the overflow and wellkilling process of a vertical H_{2}Scontaining natural gas well. The mass and momentum equations of the coupled model were solved using finite difference method, while the transient temperature prediction model were solved using finite volume method. The continuity and momentum equations are decoupled from temperature prediction model. The effect of H_{2}S content, mud displacement, drilling fluid density, and initial overflow volume on the dynamical wellkilling process of an H_{2}Scontaining natural gas well were obtained and analyzed in this work.
2 Mathematical model
2.1 Governing equations for mass and momentum
During drilling, gas enters the wellbore from the stratum when overflow occurs, and gas–liquid twophase flow appears below the location the gas reaches. Given that twophase flow exists in overflow and wellkilling periods, a mathematical model is established in this work to obtain the twophase flow behavior.
To simplify the calculation, the following assumptions are made in this work:
The gas and drilling fluid flow in a vertical wellbore is considered one dimensional;
The compressibility of the drilling fluid is ignored;
Gas and liquid phases are continuous in the control unit;
The influence of annulus eccentricity is disregarded;
The geothermal gradient is regarded as constant.
Then, the continuity equations of gas and liquid phases are expressed as follows (Wei et al., 2018):(1)(2)
On the basis of the law of conservation of momentum, the momentum equation can be obtained as follows:(3)
2.2 Transient temperature prediction model
To accurately predict the phase change behavior of H_{2}S in wellbore, a transient temperature prediction model is developed (Mao and Zhang, 2018). The unsteady twodimensional convection–diffusion and unsteady twodimensional diffusion equations are used to describe the heat transfer models as follows:
(1) Heat transfer model in the drilling string:(4)
(2) Heat transfer model of drilling string:(5)
(3) Heat transfer model in the annular:(6)
(4) Heat transfer model of casing, cement sheath, and formation:(7)
2.2 Initial and boundary conditions
(1) Overflow period
When overflow occurs, the wellhead is open and the parameters under normal drilling conditions can be used as the initial conditions of a well kick:(8)
The boundary conditions are shown as follows:(9)
(2) Killing period
After the shutin operation is performed, the flow pattern, gas void fraction, and annular pressure distributions at the shutin time can be considered the initial conditions of the wellkilling period. During the wellkilling process, the bottom hole pressure is equal to the formation pressure, and the boundary conditions are shown as follows:(10)
2.3 Flow pattern discriminant
The main patterns of a gas–liquid twophase flow are bubble, slug, churn, and annular flows (Hasan and Kabir, 1988). They can be discriminated by the following equations:
2.4 Calculation of key parameters
2.4.1 H_{2}S solubility
The solubility of H_{2}S in the drilling fluid can be expressed as (Sun et al., 2018):(16)(17)
With , , and , equation (13) can be written as(18)(19)(20)(21)(22)and(23)
Thus, the fugacity coefficient of a certain component is expressed as:(24)
The solubility of H_{2}S in the liquid can be obtained by combining equation (16) with equation (24).
2.4.2 Drift flux model
In the drift flux model, the corresponding slip rate of the gas phases and distribution coefficient in case of different flow patterns are shown as follows (Wei et al., 2018):
Thus, the gas void fraction can be calculated by the following equations:(30)(31)
2.4.3 Friction pressure drop
For single phase flow, the friction pressure drop can be obtained using the friction factor of the powerlaw fluid.(32)(33)
For gas–liquid twophase flow, the friction pressure drop can be calculated using the following equations established by Yin et al. (2017).
Churn and annular flows:(38)(39)
2.4.4 Gas production
Under the effect of pressure differential between bottom hole and formation pressures, the H_{2}Scontaining natural gas enters the wellbore. In this work, the following equation is adopted to calculate gas production (Elsharkawy, 2004):(40)(41)
2.4.5 Gas density
The gas density can be described by introducing a compression factor into the ideal gas equation, as defined in the equation below (Yin et al., 2017):(42)
2.4.6 Convective heat transfer coefficients
Convective heat transfer coefficients can be expressed as follows (Mao and Zhang, 2018):(43)
For laminar flow, Nu = 4.36; for turbulent flow, Nu = 0.023 Re^{0.8} Pr^{0.33}.
3 Model solution
3.1 Solution for mass and momentum governing equations
The finite difference method is used in this work to solve the mass and momentum governing equations. Meanwhile, the influence of the boundary and initial conditions in each control region are considered in the solving model. The solution is completed by three steps: generating discrete grids, constructing discrete equations, and solving these equations.
(1) Generating discrete grids
In finite difference numerical calculation, the spatial domain is the wellbore. For the overflow simulation, the time domain is the entire time from the overflow to the shutin. For the wellkilling simulation, the time domain is from the shutin to the time the lower boundary of twophase reaches wellhead. The schematic of discrete grid nodes of the wellbore is shown in Figure 1a.
Fig. 1 Schematic diagram of solution for the mass and momentum equations: (a) Schematic of discrete grid nodes of the wellbore; (b) Diagram of cell grid integration area K. 
(2) Model discretization
Figure 1b shows the cell grid integration area K. The partial differential equation in the mathematical model can be written as(44)
Then, by integrating this equation into area K, a curvilinear integral along the boundary of L can be obtained in accordance with Green’s theorem:(45)
The abovementioned equation can be converted to the following equation by simplification:(46)
(3) Numeralization of the continuity equation
For gasphase continuity equation, we let(47)
Thus, the gasphase difference equation can be obtained by combining equations (46) and (47):(48)
Similarly, the liquidphase difference equation can be obtained as follows:(49)
(4) Numeralization of the momentum equation
For mixed momentum equation, we let(50)
Thus, the mixed momentum difference equation can be obtained by combining equations (46) and (50):(51)
3.2 Solution for transient temperature prediction model
Figure 2 shows the grid diagram of the wellbore and formation. By using finite volume method, equations (4)–(7) can be expressed as (Mao and Zhang, 2018):(52)(53)(54)(55)
Fig. 2 The grid diagram of the wellbore and formation. 
By using underrelaxation iteration method, discretized scheme of the heat transfer control equations are obtained as follows:(56)(57)(58)(59)
Figure 3 shows the flow chart of model solution.
Fig. 3 Flow chart of the model solution process. 
3.3 Model coupling
Moreover, the continuity and momentum equations are decoupled from temperature prediction model. Consistent time step can be expressed as the following equation:(60)
When H_{2}S is released from the drilling fluid, the state changes from supercritical to gas, which directly leads to a significant increase in gas void fraction. The gasification volume of H_{2}S is influenced greatly by the temperature at a certain well depth. And the temperature also has a significant effect on the fluid properties (Xu et al., 2019). Therefore, before solving the mass and momentum equations, we must solve the transient temperature prediction model.
The wellkilling simulation is based on the overflow condition at shutin time; thus, the entire solution process can be divided into two steps: overflow and wellkilling periods. For overflow period, the solution flowchart of overflow simulation is depicted in Figure 4. For wellkilling period, the solution flowchart of wellkilling simulation is depicted in Figure 5. As shown in the figure, depending on whether the gas reaches the wellhead, the solution process can be divided into two parts: gas rising and gas exhausting stages. The boundary positions of twophase segment are key parameters in distinguishing the killing stage. The twophase flow module used in the overflow period is also called to solve the twophase flow behavior in the wellkilling period.
Fig. 4 Solution flowchart of overflow simulation for coupling model. 
Fig. 5 Solution flowchart of wellkilling simulation. 
4 Model validation
The basic data of an H_{2}Scontaining gas well Tiandong #5 were used to simulate the overflow process for verifying the validity of the dynamical wellkilling model considering an H_{2}S solubility. The overflow was found in Well Tiandong #5 at the drilling depth of 3570 m, and the pit gain was 1.65 m^{3} at 22:14. At 22:44, the wellhead emitted a loud noise, and a violent blowout accident occurred. The pit gain was recorded every 5 min, and the gas well was shut in when the pit gain reached 6 m^{3} at 22:44. Table 1 shows the basic parameters of Well Tiandong #5 (Blowout accident of TianDong #5 well, 1990). As shown in Figure 6a, the simulation results are in good agreement with the field data of Well Tiandong #5. When H_{2}S is released, the pit gain shows a sudden increase in simulation results and field data. The slight difference between the two may be due to that the compressibility of the drilling fluid is disregarded.
Fig. 6 Comparison results: (a) is the comparison between simulation and field data of Well Tiandong #5; (b) is the comparison between simulation and experimental results. 
Basic parameters of Well Tiandong #5.
Rader conducted a gas kick and wellkilling experiment in an artificial well L.S.U. 7 and provided the experimental data (Rader et al., 1975). The basic data of experimental well L.S.U. 7 were used to simulate the killing process, as shown in Table 2, for validating the validity of the wellkilling model. In the experiment, (1) nitrogen was injected into the wellbottom hole from the injection pipe until the pit gain was 1.59 m^{3}, and gas was continued to inject during the shutin operation; (2) injection of gas was stopped when the casing pressure reached the desired value; (3) well killing was started. Figure 6b shows the comparison between simulation and experimental results. As shown in Figure 6b, the simulation results of our model are in good agreement with the experimental results. The comparison between simulation results of our model and traditional “Gas column model” is also shown in Figure 6b. Large error exists in the simulation results of “Gas column model” because the friction and slip between gas and liquid are ignored.
Basic data of experimental well L.S.U 7.
5 Case study
The basic parameters from Table 3 are used to simulate the killing process of an H_{2}Scontaining natural gas well in Sichuan. The number of spatial nodes is 60 and number of spatial nodes is 80 in the simulation. The wellbore configuration of this well is shown in Figure 7. When overflow behavior occurs, the invaded gas will push the drilling fluid with the same volume out of the wellbore, thereby leading to an increase in the pit gain. In drilling engineering, once the pit gain exceeds a threshold value, a shutin operation must be performed. On the basis of field experience, the threshold value of pit gain is set as 3 m^{3} in this work.
Fig. 7 Wellbore configuration of the well in Sichuan. 
Major computational parameters of the gas well in Sichuan.
Figure 8a shows the change in gas void fraction and solubility of H_{2}S with the variation in well depth at shutin time. As shown in the figure, gas rises to a location that is 1800 m below the wellhead, and the maximal gas void fraction is approximately 0.471. The flow pattern criterion proposed by Hasan and Kabir (1988) indicates that three flow patterns exist in the wellbore: bubble, slug, and single phase flows. Besides, the H_{2}S solubility decreases as the well depth decreases. This phenomenon can be attributed to the decrease in annular temperature and pressure with the well depth, as shown in Figures 8b and 8c. When the solubility of H_{2}S decreases to H_{2}S concentration in the drilling fluid at about 500 m, the phase change of H_{2}S occurs.
Fig. 8 Change in gas void fraction, solubility of H_{2}S, annular temperature and annular pressure with the variation in well depth at shutin time: (a) Gas void fraction and solubility of H_{2}S; (b) Annular temperature; (c) Annular pressure. 
Figure 8b shows the drilling fluid temperature distribution in the drilling string and annulus at shutin time. The annular temperature is always higher than that in the drilling string because the drilling fluid in the annulus is closer to the formation, thus, more heat can be gotten from the formation. Moreover, the annular temperature and geothermal temperature intersect at about 2250 m, which indicates that when the well depth is less than 2250 m, the drilling fluid in the annulus can be cooled by the formation. When the well depth is greater than 2250 m, the drilling fluid in the annulus can be heated by the formation. Figure 8c shows the change in annular pressure with the variation in well depth at shutin time. As shown in Figure 8c, the annular pressure below 2000 m is lower than the formation pressure. This phenomenon can be attributed to an invasion of natural gas, which increases the differential pressure between bottom hole and formation. The increasing differential pressure will cause an increase in gas flow from formation to the wellbore and result in blowout accidents easily if not controlled.
Driller’s method is widely used in well killing in oil and gas fields (Feng et al., 2016a). Thus, the method is also used in the simulation of this work. Figure 9 shows the physical model of the killing process. During killing, the pressure balance relationship can be expressed as follows (Feng et al., 2016b):(61)
Fig. 9 Physical model of the killing process. 
Figure 10a shows the change in bottom hole pressure with the variation in time. As shown in the figure, the bottom hole pressure first decreases during the overflow period and then keeps constant (greater or equal to formation pressure) when well killing is conducted. Therefore, no gas flows into the wellbore from the formation during killing.
Fig. 10 Change in bottom hole pressure and casing pressure the variation in time: (a) Bottom hole pressure; (b) Casing pressure. 
With the change in hydrostatic pressure created by the fluid column, the casing pressure of the wellhead can be changed by adjusting the throttle opening to maintain the bottom hole pressure constant. Figure 10b shows the change in casing pressure with the variation in time. As shown in the figure, the entire killing process can be divided into six stages. Stage 1 is the initial stage of the killing process, which is determined by the flow pattern at shutin time. Stage 2 indicates that the twophase section in the wellbore is rising up toward the wellbore prior to reaching the wellhead. At this time, the flow pattern distribution in the wellbore is single phase flow, gas–liquid twophase flow, and single phase flow in sequence. Moreover, the casing pressure increases with the increase in killing time. This phenomenon is attributed to the slippage and expansion of gas. The gas will expand as it rises up along the wellbore because the annular and formation temperatures and depth decrease (Sun et al., 2018). Thus, the hydrostatic pressure decreases accordingly and the casing pressure increases to keep the bottom hole pressure constant. Stage 3 indicates that the gas reaches the wellhead and single phase and gas–liquid twophase flows exist in the wellbore. The total gas volume reaches a maximum at this point. Stage 4 indicates that all the gas has been pushed out of the wellbore and only single phase flow exists in the wellbore. This stage will last for some time until the weighted drilling fluid enters the annulus. Furthermore, the casing pressure keeps constant in Stage 4. Stage 5 indicates that the weighted drilling fluid enters the annulus and the hydrostatic pressure increases. Thus, the casing pressure decreases. Stage 6 indicates that the annulus is full of the weighted drilling fluid and the casing pressure decreases to 0. The wellkilling operation ends.
Figure 10b also shows that, when all the gas is about to be pushed out of the well, H_{2}S is released, and the curve drops slowly. This phenomenon is the result of casing pressure applied at the wellhead. Figure 8a shows that, when overflow behavior occurs, H_{2}S will release at a location close to the wellhead if shutin operation is not performed. However, H_{2}S solubility decreases with the increase in pressure, and the casing pressure applied at the wellhead changes the pressure distribution in the wellbore (Fig. 11). As shown in Figure 11, the pressure of each location along the wellbore is higher than that of the overflow process during the killing process. Therefore, the H_{2}S solubility of each location in the killing process increases compared with that of shutin time. Although it reaches the wellhead, the H_{2}S will not release until the H_{2}S solubility is lower than its content in the drilling fluid. Therefore, when the casing pressure decreases, the H_{2}S solubility decreases. As a result, the H_{2}S transforms from a supercritical phase to a gas phase, gas volume expands suddenly, and gas void fraction increases abruptly. Thus, the casing pressure drops slowly.
Fig. 11 Dynamical pressure distribution in the wellbore: (a) during gas kick; (b) during well killing. 
5.1 Analysis under different H_{2}S contents
Considering that the phase change of H_{2}S occurs at low pressure and temperature and the phase change of H_{2}S has a great effect on the pressure distribution (Sun et al., 2013), the diagrams of wellhead flow pattern when gas reaches wellhead during well killing under different working conditions are depicted to study the flows at wellhead.
Figure 12 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with H_{2}S contents of 0%, 2%, 4%, and 6%. As shown in the figure, the gas void fraction increases with the increase in H_{2}S content. This phenomenon is caused by the increase in gasification volume of H_{2}S as H_{2}S content increases. When the well depth decreases, the annular pressure and annular temperature decrease, as shown in Figure 13. Besides, the higher the H_{2}S content is, the lower the annular temperature will be as shown in Figure 13b. This is because, near the critical point of H_{2}S, the viscosity of supercritical H_{2}S decreases significantly (Vesovic et al., 1998). Therefore, in the upper well section, the convective heat transfer coefficient between the drilling fluid and the formation increases as viscosity of supercritical H_{2}S decreases. Thus, more heat will be diffused from the drilling fluid in the annulus to the formation. As H_{2}S content increases, this effect will be enhanced. Therefore, in the upper well section which is near the critical point of H_{2}S, the annular temperature is lower at higher H_{2}S content and, as a result, at the same well depth, the H_{2}S solubility will be lower and more H_{2}S will gasify.
Fig. 12 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with H_{2}S contents of 0%, 2%, 4%, and 6%. 
Fig. 13 Change in annular pressure and temperature with the variation in well depth at 1300 s during well killing with H_{2}S contents of 0%, 2%, 4%, and 6%: (a) Annular pressure; (b) Annular temperature. 
Additional H_{2}S will gasify at higher H_{2}S content, and this condition changes the wellhead flow pattern. Churn flow appears at the wellhead with H_{2}S content of 0%; however, when H_{2}S is released, annular flow appears and the droplets in the bubble decrease as H_{2}S content increases. Therefore, the hydrostatic pressure decreases with the increase in H_{2}S content, and the casing pressure increases to keep the bottom hole pressure constant. Figure 14 shows that, when H_{2}S is released, the casing pressure increases with H_{2}S content. Thus, the annular pressure of the upper well section increases with the increase in H_{2}S content, as shown in Figure 13a. However, H_{2}S content has a small effect on entire wellbore pressure distribution, as shown in Figure 14.
Fig. 14 Change in wellbore pressure and casing pressure with the variation in time and H_{2}S contents of 0%, 2%, 4%, and 6%: 
From comparison between Figures 8c and 13a, it can be concluded that well killing can rebuild the pressure distribution in the wellbore and stop the overflow behavior caused by the pressure differential between bottom hole and formation. However, for the stratum with low formation fracture pressure, the wellkilling operation may cause annulus pressure to approach or exceed formation fracture pressure, thereby causing well leakage. This scenario should be highly regarded, especially for H_{2}Scontaining natural gas wells. Figure 13a shows that, when H_{2}S is released, the annular pressure of the upper well section will exceed the formation fracture pressure. Moreover, a high H_{2}S content indicates serious fracture of the stratum. Therefore, during the late period of well killing in H_{2}Scontaining natural gas wells, the casing pressure should be reduced appropriately to avoid stratum fracturing and prevent leakage accidents.
5.2 Analysis under different initial overflow volumes
Figure 15 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with initial overflow volumes of 2, 4, and 8 m^{3}. As shown in the figure, the gas void fraction increases with the increase in initial overflow volume because a large amount of gas enters the wellbore at high initial overflow volume. Moreover, the wellhead flow pattern with initial overflow volume of 2 m^{3} can be regarded as slug flow, whereas the annular flow appears at high initial overflow volumes of 4 and 8 m^{3}. The gas phase becomes continuous at initial overflow volumes of 4 and 8 m^{3}. The continuous phase volume increases as the initial overflow volume increases.
Fig. 15 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with initial overflow volumes of 2, 4, and 8 m^{3}. 
A high initial overflow volume indicates gas rises to a high location in the wellbore at shutin time. When initial overflow volume reaches 8 m^{3}, the gas enters the wellhead at shutin time. This conclusion can also be drawn from Figure 16. Figure 16b shows the change in casing pressure with the variation in time and initial overflow volumes of 2, 4, and 8 m^{3}. As shown in the figure, the surface casing pressure increases with the increase in initial overflow volume. This phenomenon can be attributed to the increase in hydrostatic pressure loss as gas volume increases. Thus, additional surface casing pressure should be applied to balance the bottom hole pressure at high initial overflow volume. Moreover, as shown in Figure 16b, the pressure change in the wellbore intensifies as initial overflow volume increases.
Fig. 16 Change in wellbore pressure and casing pressure with the variation in time and initial overflow volumes of 2, 4, and 8 m^{3}: (a) Wellbore pressure; (b) Casing pressure. 
Besides, as the initial overflow volume increases, the gasification starting time of H_{2}S increases and, the casing pressure is lower as shown in Figure 16a. This phenomenon can be attributed to the change of drilling fluid temperature. Figure 17a shows the change in annular temperature with the variation in well depth at 1300 s during well killing with different initial overflow volume. We can see that, the drilling fluid temperature is higher at higher initial overflow volume. This is because, as more gas enters the wellbore, more heat from the formation will be carried into the wellbore. Thus, the H_{2}S solubility will be higher at the same well depth and less H_{2}S will gasify. Therefore, when H_{2}S gasifies, the casing pressure applied at wellhead is higher at lower initial overflow volume.
Fig. 17 Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with initial overflow volumes of 2, 4, and 8 m^{3}: (a) Annular temperature; (b) Annular pressure. 
Figure 17b shows the change in annular pressure with the variation in well depth and initial overflow volumes of 2, 4, and 8 m^{3}. As initial overflow volume increases, the gas expansion capacity is limited. Thus, pressure changed rapidly from 2 m^{3} to 4 m^{3}, but slowly from 4 m^{3} to 8 m^{3}. As shown in the figure, the annular pressure is high with high initial overflow volume. This phenomenon is caused by the increase in surface casing pressure as initial overflow volume increases. However, when initial overflow volume increases, the annular pressure may approach or even exceed the formation fracture pressure, which may lead to formation fracturing. Therefore, careful monitoring of the pit gain is necessary during the exploitation of H_{2}Scontaining natural gas wells. Once the pit gain exceeds a threshold value, the shutin operation should be performed immediately to prevent further overflow. Furthermore, during the late period of well killing in H_{2}Scontaining natural gas wells with high initial overflow volume, the surface casing pressure should be decreased appropriately to avoid stratum fracturing and prevent leakage accidents. Early detection of gas kick should be more frequent to avoid severe overflow.
5.3 Analysis under different drilling fluid densities
Figure 18 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}. As shown in the figure, the gas void fraction at wellhead increases with the increase in drilling fluid density. This phenomenon can be attributed to the decrease in casing pressure caused by the increase in drilling fluid density. The annular pressure is low at shutin time with low drilling fluid density because the hydrostatic pressure is proportional to drilling fluid density. Therefore, to maintain the bottom hole pressure, the surface casing pressure should be higher at low drilling fluid density, as shown in Figure 19. During well killing, the annular pressure is high at low drilling fluid density, as shown in Figure 20a. The gas void fraction is high at high drilling fluid density because the gas expands violently at low pressure. Figure 18 shows that the drilling fluid density slightly affects the wellhead flow pattern, and annular flow is the common flow pattern under drilling fluid densities of 1.20, 1.25, and 1.30 g/m^{3}.
Fig. 18 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}. 
Fig. 19 Change in wellbore pressure and casing pressure with the variation in time and drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}: (a) Wellbore pressure; (b) Casing pressure. 
Fig. 20 Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}: (a) Annular pressure; (b) Annular temperature. 
Figure 20b shows the change in annular temperature with the variation in well depth at 1300 s during well killing with different drilling fluid density. As shown in Figure 20b, for the deeper well section, the drilling fluid temperature increases with the decrease in drilling fluid density. This is because, when the same heat is obtained from the formation, the greater the drilling fluid density, the less the increase in temperature (Gao et al., 2017). For the upper well section, the drilling fluid temperature is closer to the formation temperature, the temperature differences decrease. Considering that H_{2}S only gasifies in the upper well section, the effect of temperature difference caused by the change in drilling fluid density on gasification of H_{2}S is small as shown in Figure 19b.
Furthermore, annular pressure of the upper well section will exceed the formation fracture pressure, and low drilling fluid density indicates serious fracture of the stratum. Thus, during the well killing in H_{2}Scontaining natural gas wells, higher drilling fluid density should be used to avoid stratum fracturing and prevent leakage accidents.
5.4 Analysis under different wellkilling displacements
Figure 21 shows the diagram of wellhead flow pattern when gas reaches wellhead during well killing with displacements of 20, 30, and 40 L/s. As shown in the figure, the gas void fraction at wellhead increases slightly with the increase in displacement. The annular flow is the common wellhead flow pattern under displacements of 20, 30, and 40 L/s, and the number of droplets in the continuous gas decreases with the increase in displacement.
Fig. 21 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with displacements of 20, 30, and 40 L/s. 
Figure 22a shows the change in annular temperature with the variation in well depth at 1300 s during well killing with different displacement. The lower the displacement is, the longer time for heat exchange with the formation will be. Thus, for the upper well section, the annular temperature is higher than formation, at higher displacement, the drilling fluid temperature is higher because lesser heat of drilling fluid will loss. For the deeper well section, the annular temperature is lower than formation, at higher displacement, the drilling fluid temperature is lower because lesser heat from formation will be exchanged to the drilling fluid. However, the effect of displacement on temperature distribution in annulus is small.
Fig. 22 Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with displacements of 20, 30, and 40 L/s: (a) Annular temperature; (b) Annular pressure. 
Figure 23 shows the change in wellbore pressure and casing pressure with the variation in time and displacements of 20, 30, and 40 L/s. As shown in the figure, the higher the displacement is, the faster the gas kick leaves wellbore and, the lower the casing pressure will be. Besides, H_{2}S will gasify at an earlier time. This phenomenon can be attributed to the increased friction resistance due to the increase in displacement (Shi et al., 2014). Thus, the bottom hole pressure increases as displacement increases and, the casing pressure decreases as a result. Therefore, the annular pressure is lower at higher displacement in the upper well section as shown in Figure 22b. Besides, the annular pressure of the upper well section will exceed the formation fracture pressure, and low displacement indicates serious fracture of the stratum. Thus, during well killing in H_{2}Scontaining natural gas wells, higher wellkilling displacement should be used to avoid stratum fracturing and prevent leakage accidents.
Fig. 23 Change in wellbore pressure and casing pressure with the variation in time with displacements of 20, 30, and 40 L/s: (a) Wellbore pressure; (b) Casing pressure. 
6 Conclusion
This work established a dynamical wellkilling model considering an H_{2}S solubility to simulate the wellkilling process of a vertical H_{2}Scontaining natural gas well. The following important conclusions can be drawn:
H_{2}S will gasify near wellhead during well killing when casing pressure decreases. Near the critical point of H_{2}S, the annular temperature decreases as H_{2}S content increases and, the H_{2}S solubility will be lower and more H_{2}S will gasify. To balance the bottom hole pressure, when H_{2}S releases, the casing pressure increases as H_{2}S content increase.
As initial overflow volume increases, the annular temperature, annular pressure and the casing pressure increase significantly. When H_{2}S gasifies, the casing pressure applied at wellhead should be higher at lower initial overflow volume to balance bottom hole pressure.
In the wellkilling process, the annular pressure and temperature decrease as drilling fluid density increases and a lower casing pressure is needed for balancing bottom hole pressure. No significant effect of drilling fluid density on the gasification of H_{2}S is found. As wellkilling displacement increases, for the upper well section, the annular temperature increases while the annular pressure decreases. The casing pressure is lower at a higher displacement for higher friction resistance. Besides, H_{2}S will gasify at an earlier time.
When drilling for H_{2}Scontaining natural gas well, early detection of gas kick should be more frequent to avoid severe overflow. Besides, during well killing, higher displacement and density of drilling fluid should be considered to avoid stratum fracturing and prevent leakage accidents under the premise of meeting drilling requirements.
Acknowledgments
This work was supported by the National Key Research and Development Program of China (2019YFC0312303) and Sichuan Science and Technology Project (2019YFS0045).
References
 AmayaGómez R., López J., Pineda H., UrbanoCaguasango D., Pinilla J., Ratkovich N., Muñoz F. (2019) Probabilistic approach of a flow pattern map for horizontal, vertical, and inclined pipes, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 74, 67. [CrossRef] [Google Scholar]
 Ansari A.M., Sylvester N.D., Sarica C., Shoham O., Brill J.P. (1994) A comprehensive mechanistic model for upward two phase flow in wellbores, SPE Prod. Facil. 5, 143–152. https://doi.org/10.2523/20630MS. [CrossRef] [Google Scholar]
 Bilicki Z., Kestin J. (1987) Transition criteria for twophase flow patterns in vertical upward flow, Int. J. Multiphase Flow. 13, 3, 283–294. [CrossRef] [Google Scholar]
 Blowout accident of TianDong well #5 (1990) Petrochina Southwest Oil and Gas field Company. [Google Scholar]
 Ebrahimi A., Khamehchi E. (2015) A robust model for computing pressure drop in vertical multiphase flow, J. Nat. Gas. Sci. Eng. 26, 1306–1316. [Google Scholar]
 Elsharkawy A.M. (2004) Efficient methods for calculations of compressibility, density and viscosity of natural gases, Fluid Phase Equilibr. 218, 1, 1–13. [CrossRef] [Google Scholar]
 Feng J., Fu J., Chen P., Liu Z., Wei H. (2015) Predicting pressure behavior during dynamic kill drilling with a twophase flow, J. Nat. Gas Sci. Eng. 22, 591–597. [Google Scholar]
 Feng J., Fu J., Chen P., Luo J., Kou B. (2016a) Comparisons of the Driller’s method and the wait and weight method in deepwater well killing operation, Arab. J. Sci. Eng. 41, 7, 2699–2706. [Google Scholar]
 Feng J., Fu J., Luo J., Liu Z., Wei H. (2016b) An advanced Driller’s Method simulator for deepwater well control, J. Loss Prevent. Proc. 39, 131–140. [CrossRef] [Google Scholar]
 Fu J., Su Y., Jiang W., Li S., Chen Y. (2019) Wellbore annulus water hammer pressure prediction based on transient multiphase flow characteristics, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 74, 84. [CrossRef] [Google Scholar]
 Gao Y., Cui Y., Xu B., Sun B., Zhao X., Li H., Chen L. (2017) Two phase flow heat transfer analysis at different flow patterns in the wellbore, Appl. Therm. Eng. 117, 544–552. [Google Scholar]
 Guo R., Chen Y., Waltrich P.J., Williams W.C. (2018) An experimental investigation on flow pattern map and driftflux model for CoCurrent upward liquidgas twophase flow in narrow annuli, J. Nat. Gas. Sci. Eng. 51, 65–72. [Google Scholar]
 Guo X., Wang Q. (2016) A new prediction model of elemental sulfur solubility in sour gas mixtures, J. Nat. Gas Sci. Eng. 31, 98–107. [Google Scholar]
 Hasan A.R., Kabir C.S. (1988) A study of multiphase flow behavior in vertical wells, Spede. 3, 2, 263–272. https://doi.org/10.2118/15138PA. [Google Scholar]
 Hou Z., Yan T., Li Z., Feng J., Sun S., Yuan Y. (2019) Temperature prediction of two phase flow in wellbore using modified heat transfer model: An experimental analysis, Appl Therm. Eng. 149, 54–61. [Google Scholar]
 Kelessidis V.C., Dukler A.E. (1989) Modeling flow pattern transitions for upward gasliquid flow in vertical concentric and eccentric annuli, Int. J. Multiphase Flow. 15, 2, 173–191. [CrossRef] [Google Scholar]
 Li H.Z., Mouline Y. (1997) Chaotic bubble coalescence in nonnewtonian fluids, Int. J. Multiphase Flow. 23, 713–723. [CrossRef] [Google Scholar]
 Mao L., Zhang Z. (2018) Transient temperature prediction model of horizontal wells during drilling shale gas and geothermal energy, J. Petrol. Sci. Eng. 169, 610–622. [CrossRef] [Google Scholar]
 Meng Y., Xu C., Na W., Gao L., Hongtao L., Mubai D. (2015) Numerical simulation and experiment of the annular pressure variation caused by gas kick/injection in wells, J Nat Gas Sci Eng. 22, 646–655. [Google Scholar]
 Nickens H.V. (1987) A dynamic computer model of kick Well, SPE Drill. Eng. 2, 2, 158–173. https://doi.org/10.2118/14183PA. [Google Scholar]
 Osman E.S.A. (2004) Artificial neural network models for identifying flow regimes and predicting liquid holdup in horizontal multiphase flow, SPE Prod. Facil. 19, 1, 33–40. https://doi.org/10.2118/68219MS. [CrossRef] [Google Scholar]
 Rader D.W., Bourgoyne A.T., Ward R.H. (1975, May 1) Factors affecting bubblerise velocity of gas kicks, Society of Petroleum Engineers. https://doi.org/10.2118/4647PA. [Google Scholar]
 Raghavan R. (1989) Welltest analysis for multiphase flow, SPE Form. Eval. 4, 4, 585–594. https://doi.org/10.2118/14098PA. [CrossRef] [Google Scholar]
 Shi H., Li G., Huang Z., Shi S. (2014) Properties and testing of a hydraulic pulse jet and its application in offshore drilling, Petrol. Sci. 11, 3, 401–407. [CrossRef] [Google Scholar]
 Sun B., Gong P.B., Wang Z.Y. (2013) Simulation of gas kick with high H_{2}S content in deep well, J. Hydrodyn. 25, 2, 264–273. [CrossRef] [Google Scholar]
 Sun B.J., Sun X.H., Wang Z.Y., Chen Y.H. (2017) Effects of phase transition on gas kick migration in deepwater horizontal drilling, J. Nat. Gas Sci. Eng. 46, 710–729. [Google Scholar]
 Sun B.J., Guo Y.H., Sun W.T., Gao Y., Hao L., Wang Z., Zhang H. (2018) Multiphase flow behavior for acidgas mixture and drilling fluid flow in vertical wellbore, J. Petrol. Sci. Eng. 165, 388–396. [CrossRef] [Google Scholar]
 Vesovic V., Assael M.J., Gallis Z.A. (1998) Prediction of the viscosity of supercritical fluid mixtures, Int. J. Thermophys. 19, 5, 1297–1313. [Google Scholar]
 Wang N., Sun B., Wang Z., Wang J., Yang C. (2016) Numerical simulation of two phase flow in wellbores by means of drift flux model and pressure based method, J. Nat. Gas. Sci. Eng. 36, 811–823. [Google Scholar]
 Wang Z.Y., Sun B.J. (2014) Deepwater gas kick simulation with consideration of the gas hydrate phase transition, J. Hydrodyn. 26, 1, 94–103. [CrossRef] [Google Scholar]
 Wang Z., Sun B. (2009) Annular multiphase flow behavior during deep water drilling and the effect of hydrate phase transition, Petrol. Sci. 6, 1, 57–63. [CrossRef] [Google Scholar]
 Wei N., Xu C.Y., Meng Y.F., Li G. (2018) Numerical simulation of gasliquid twophase flow in wellbore based on drift flux model, Appl. Math. Comput. 338, 175–191. [Google Scholar]
 Xu Z.M., Song X.Z., Li G.S., Wu K., Pang Z., Zhu Z. (2018) Development of a transient nonisothermal twophase flow model for gas kick simulation in HTHP deep well drilling, Appl. Therm. Eng. 141, 1055–1069. [Google Scholar]
 Xu Z.M., Song X.Z., Li G.S., Zhu Z., Zhu B. (2019) Gas kick simulation in oilbased drilling fluids with the gas solubility effect during hightemperature and highpressure well drilling, Appl. Therm. Eng. 149, 1080–1097. [Google Scholar]
 Xie J., Zhang X., Tang Y., Wang Y., Shao Q., Yu B. (2014) Transient simulation of the blowingout process of the air pockets in vertical wellbore, Appl. Therm. Eng. 72, 97–103. [Google Scholar]
 Yin H., Liu P., Li Q., Wang Q., Gao D. (2015) A new approach to risk control of gas kick in highpressure sour gas wells, J. Nat. Gas Sci. Eng. 26, 142–148. [Google Scholar]
 Yin B., Gang L., Li X. (2017) Multiphase transient flow model in wellbore annuli during gas kick in deepwater drilling based on oilbased mud, Appl. Math. Model. 51, 159–198. [Google Scholar]
 Zhang X., Li Q., Zheng L., Li X., Xu L. (2020) Numerical simulation and feasibility assessment of acid gas injection in a carbonate formation of the Tarim Basin, China. Oil Gas, Sci. Technol. 75, 1–17. [Google Scholar]
 Zhang J., Li X., Tang X., Luo W. (2019) Establishment and analysis of temperature field of riserless mud recovery system, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 74, 19, 1–9. [CrossRef] [Google Scholar]
 Zhu G., Zhang S., Huang H., Liang Y., Meng S., Li Y. (2011) Gas genetic type and origin of hydrogen sulfide in the Zhongba gas field of the western Sichuan Basin, China, Appl. Geochem. 26, 7, 1261–1273. [Google Scholar]
All Tables
All Figures
Fig. 1 Schematic diagram of solution for the mass and momentum equations: (a) Schematic of discrete grid nodes of the wellbore; (b) Diagram of cell grid integration area K. 

In the text 
Fig. 2 The grid diagram of the wellbore and formation. 

In the text 
Fig. 3 Flow chart of the model solution process. 

In the text 
Fig. 4 Solution flowchart of overflow simulation for coupling model. 

In the text 
Fig. 5 Solution flowchart of wellkilling simulation. 

In the text 
Fig. 6 Comparison results: (a) is the comparison between simulation and field data of Well Tiandong #5; (b) is the comparison between simulation and experimental results. 

In the text 
Fig. 7 Wellbore configuration of the well in Sichuan. 

In the text 
Fig. 8 Change in gas void fraction, solubility of H_{2}S, annular temperature and annular pressure with the variation in well depth at shutin time: (a) Gas void fraction and solubility of H_{2}S; (b) Annular temperature; (c) Annular pressure. 

In the text 
Fig. 9 Physical model of the killing process. 

In the text 
Fig. 10 Change in bottom hole pressure and casing pressure the variation in time: (a) Bottom hole pressure; (b) Casing pressure. 

In the text 
Fig. 11 Dynamical pressure distribution in the wellbore: (a) during gas kick; (b) during well killing. 

In the text 
Fig. 12 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with H_{2}S contents of 0%, 2%, 4%, and 6%. 

In the text 
Fig. 13 Change in annular pressure and temperature with the variation in well depth at 1300 s during well killing with H_{2}S contents of 0%, 2%, 4%, and 6%: (a) Annular pressure; (b) Annular temperature. 

In the text 
Fig. 14 Change in wellbore pressure and casing pressure with the variation in time and H_{2}S contents of 0%, 2%, 4%, and 6%: 

In the text 
Fig. 15 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with initial overflow volumes of 2, 4, and 8 m^{3}. 

In the text 
Fig. 16 Change in wellbore pressure and casing pressure with the variation in time and initial overflow volumes of 2, 4, and 8 m^{3}: (a) Wellbore pressure; (b) Casing pressure. 

In the text 
Fig. 17 Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with initial overflow volumes of 2, 4, and 8 m^{3}: (a) Annular temperature; (b) Annular pressure. 

In the text 
Fig. 18 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}. 

In the text 
Fig. 19 Change in wellbore pressure and casing pressure with the variation in time and drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}: (a) Wellbore pressure; (b) Casing pressure. 

In the text 
Fig. 20 Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with drilling fluid densities of 1.20, 1.25, and 1.30 g/cm^{3}: (a) Annular pressure; (b) Annular temperature. 

In the text 
Fig. 21 Diagram of wellhead flow pattern when gas reaches wellhead during well killing with displacements of 20, 30, and 40 L/s. 

In the text 
Fig. 22 Change in annular temperature and annular pressure with the variation in well depth at 1300 s during well killing with displacements of 20, 30, and 40 L/s: (a) Annular temperature; (b) Annular pressure. 

In the text 
Fig. 23 Change in wellbore pressure and casing pressure with the variation in time with displacements of 20, 30, and 40 L/s: (a) Wellbore pressure; (b) Casing pressure. 

In the text 