Dynamical well-killing simulation of a vertical H2S-containing natural gas well

This work aims to explore the dynamical well-killing process of a vertical H2S-containing natural gas well. A dynamical well-killing model considering an H2S solubility was established to simulate the overflow and well-killing process of a vertical H2S-containing 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 H2S content, mud displacement, drilling fluid density, and initial overflow volume on the dynamical well-killing process of an H2S-containing natural gas well were obtained and analyzed in this work. Results showed that H2S will gasify near wellhead during well killing when casing pressure decreases. To balance the bottom hole pressure, when H2S releases, the casing pressure increases as H2S content increases. As initial overflow volume increases, the annular temperature, annular pressure and the casing pressure increase significantly. When H2S gasifies, the casing pressure applied at wellhead should be higher at lower initial overflow volume to balance bottom hole pressure. In the well-killing 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 well-killing displacement increases H2S will gasify at an earlier time. When drilling for H2S-containing 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.


Q g
Gas production Q gÀH 2 S Released H 2 S q g Densities of gas q l Density of drilling fluid v g Velocity of gas v l Velocity of drilling fluid E g Gas void fraction E l Liquid holdup A Cross-sectional area g Local acceleration of gravity

Introduction
Gas kick and well-killing 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 (Amaya-Gó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 two-phase 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 . Their model has also been used to study the phase transition of hydrate during the exploitation of natural gas hydrate Sun, 2009, 2014). Xu et al. (2018Xu et al. ( , 2019) developed a non-isothermal two-phase flow model to investigate the effect of major parameters on the two-phase 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 gas-liquid two-phase 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 multi-phase 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 two-phase flow models may be unsuitable for predicting the multiphase flow behavior in H 2 S-containing 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 S-containing 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(Sun et al., , 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 S-containing natural gas well. By applying casing pressure in the wellhead during well killing, the two-phase 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 well-killing process of an H 2 S-containing natural gas well to offer a guidance for exploitation of H 2 S-containing gas reservoir.
Thus, the aim of this study is to establish a dynamical well-killing model considering an H 2 S solubility to simulate the overflow and well-killing 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 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 well-killing process of an H 2 S-containing natural gas well were obtained and analyzed in this work.

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 two-phase flow appears below the location the gas reaches. Given that two-phase flow exists in overflow and well-killing periods, a mathematical model is established in this work to obtain the two-phase flow behavior.
To simplify the calculation, the following assumptions are made in this work: 1. The gas and drilling fluid flow in a vertical wellbore is considered one dimensional; 2. The compressibility of the drilling fluid is ignored; 3. Gas and liquid phases are continuous in the control unit; 4. The influence of annulus eccentricity is disregarded; 5. The geothermal gradient is regarded as constant.
Then, the continuity equations of gas and liquid phases are expressed as follows (Wei et al., 2018): On the basis of the law of conservation of momentum, the momentum equation can be obtained as follows:

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: (2) Heat transfer model of drilling string: (3) Heat transfer model in the annular: (4) Heat transfer model of casing, cement sheath, and formation:

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: The boundary conditions are shown as follows: (2) Killing period After the shut-in operation is performed, the flow pattern, gas void fraction, and annular pressure distributions at the shut-in time can be considered the initial conditions of the well-killing period. During the well-killing process, the bottom hole pressure is equal to the formation pressure, and the boundary conditions are shown as follows:

Flow pattern discriminant
The main patterns of a gas-liquid two-phase flow are bubble, slug, churn, and annular flows (Hasan and Kabir, 1988). They can be discriminated by the following equations: Bubble flow: Slug flow: Churn flow: v sg < 3:1 Annular flow: v sg > 3:1 gðq l À q g Þr q 2 g " # 0:25 : 2.4 Calculation of key parameters

H 2 S solubility
The solubility of H 2 S in the drilling fluid can be expressed as (Sun et al., 2018): With A ¼ ap R 2 T 2 , B ¼ bp RT , and Z ¼ PV RT , equation (13) can be written as and m ¼ 0:37464 þ 1:5226 À 0:26992w 2 : Thus, the fugacity coefficient of a certain component is expressed as: The solubility of H 2 S in the liquid can be obtained by combining equation (16) with equation (24).

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): Bubble flow: v gr ¼ 1:53 gðq l À q g Þr q l 2 0:25 ; ð25Þ Slug flow: v gr ¼ 0:35 gDðq l À q g Þ q l 0:5 ; ð27Þ Churn and annular flows: Thus, the gas void fraction can be calculated by the following equations:

Friction pressure drop
For single phase flow, the friction pressure drop can be obtained using the friction factor of the power-law fluid.
For gas-liquid two-phase flow, the friction pressure drop can be calculated using the following equations established by Yin et al. (2017).

Gas production
Under the effect of pressure differential between bottom hole and formation pressures, the H 2 S-containing natural gas enters the wellbore. In this work, the following equation is adopted to calculate gas production (Elsharkawy, 2004):

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): 2.4.6 Convective heat transfer coefficients Convective heat transfer coefficients can be expressed as follows (Mao and Zhang, 2018): For laminar flow, Nu = 4.36; for turbulent flow, Nu = 0.023 Re 0.8 Pr 0.33 .

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 shut-in. For the well-killing simulation, the time domain is from the shut-in to the time the lower boundary of two-phase reaches wellhead. The schematic of discrete grid nodes of the wellbore is shown in Figure 1a.
(2) Model discretization Figure 1b shows the cell grid integration area K. The partial differential equation in the mathematical model can be written as 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: The above-mentioned equation can be converted to the following equation by simplification: (3) Numeralization of the continuity equation For gas-phase continuity equation, we let Thus, the gas-phase difference equation can be obtained by combining equations (46) and (47): Similarly, the liquid-phase difference equation can be obtained as follows: (4) Numeralization of the momentum equation For mixed momentum equation, we let Thus, the mixed momentum difference equation can be obtained by combining equations (46) and (50): 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): By using under-relaxation iteration method, discretized scheme of the heat transfer control equations are obtained as follows: See equations (56)-(59) top of next page Figure 3 shows the flow chart of model solution.

Model coupling
Moreover, the continuity and momentum equations are decoupled from temperature prediction model. Consistent time step can be expressed as the following equation: 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 Schematic of discrete grid nodes of the wellbore. 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 well-killing simulation is based on the overflow condition at shut-in time; thus, the entire solution process can be divided into two steps: overflow and well-killing periods. For overflow period, the solution flowchart of overflow simulation is depicted in Figure 4. For well-killing period, the solution flowchart of well-killing 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 two-phase segment are key parameters in distinguishing the killing stage. The twophase flow module used in the overflow period is also called to solve the two-phase flow behavior in the well-killing period.

Model validation
The basic data of an H 2 S-containing gas well Tiandong #5 were used to simulate the overflow process for verifying the validity of the dynamical well-killing 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. Rader conducted a gas kick and well-killing 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 well bottom hole from the injection pipe until the pit gain was 1.59 m 3 , and gas was continued to inject during the shut-in 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.

Case study
The basic parameters from Table 3 are used to simulate the killing process of an H 2 S-containing 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 shut-in 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.    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. Figure 8b shows the drilling fluid temperature distribution in the drilling string and annulus at shut-in 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 shut-in 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): 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.
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 shut-in time. Stage 2 indicates that the two-phase 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 two-phase 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 two-phase 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 well-killing 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 shut-in 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 shut-in 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.

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.
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.
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 well-killing 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 S-containing 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 S-containing natural gas wells, the casing pressure should be reduced appropriately to avoid stratum fracturing and prevent leakage accidents. 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.

Analysis under different initial overflow volumes
A high initial overflow volume indicates gas rises to a high location in the wellbore at shut-in time. When initial overflow volume reaches 8 m 3 , the gas enters the wellhead at shut-in 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.    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. 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 S-containing natural gas wells. Once the pit gain exceeds a threshold value, the shut-in operation should be performed immediately to prevent further overflow. Furthermore, during the late period of well killing in H 2 S-containing 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. 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 shut-in 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 . 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.

Analysis under different drilling fluid densities
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 S-containing natural gas wells, higher drilling fluid density should be used to avoid stratum fracturing and prevent leakage accidents. 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. 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.

Analysis under different well-killing displacements
However, the effect of displacement on temperature distribution in annulus is small. 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 S-containing natural gas wells, higher wellkilling displacement should be used to avoid stratum fracturing and prevent leakage accidents.

Conclusion
This work established a dynamical well-killing model considering an H 2 S solubility to simulate the well-killing process of a vertical H 2 S-containing natural gas well. The following important conclusions can be drawn: 1. 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. 2. 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. 3. In the well-killing 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 well-killing 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. 4. When drilling for H 2 S-containing 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.