A Multiphase Wellbore Flow Model for Sour Gas “ Kicks ”

: This study presents a new multiphase ﬂ ow model with transient heat transfer and pressure coupling to simulate HTHP (high temperature and high pressure) sour gas “ kicks ” phenomena. The model is intended to support the estima-tion of wellbore temperature and pressure when sour gas kicks occur during drilling operation. The model considers sour gas solubility, phase transition and effects of temperature and pressure on the physical parameters of drilling ﬂ uid. Experimental data for a large-diameter pipe ﬂ ow are used to validate the model. The results indicate that with ﬂ uid circulation, the annulus temperature with H 2 S kicks is the highest, followed by CO 2 , and CH 4 is the lowest. The phase transition point of H 2 S is closer to wellhead compared with CO 2 , resulting in a faster expansion rate, which is more imperceptible and dangerous. With ﬂ uid circulation, the drilling ﬂ uid density and plastic viscosity both ﬁ rst decrease and then increase with the increase in the well depth. The bottom hole pressure when H 2 S kicks is greater than that for CO 2 with the same amount of sour gas, and the pressure difference gradually increases with the increase of H 2 S/CO 2 content. In addition, a parametric sensitivity analysis has been conducted to evaluate qualitatively and rank the in ﬂ uential factors affecting the bottom hole temperature and pressure.


Introduction
Widely distributed in Sichuan Basin and Tarim Basin in China, as well as in some other areas around the world, sour gas reservoirs are becoming one of the key areas in oil and gas exploration and development. About 2/3 of the gas fields in Sichuan Basin in China contain hydrogen sulfide [1]. The highest CO 2 content in the center of Tarim Basin in China reaches 67% [2], and the H 2 S content of gas fields at South Texas in the U.S. is up to 98%. Sour gas reservoirs are generally located in deep marine carbonate formation and are characterized by high temperature and high pressure (HTHP), complex formation pressure system and narrow mud density window, which frequently cause complicated accidents during drilling operations [3]. Sour gas could easily invade into wellbore and cause kicks and blowouts while drilling sour gas reservoirs especially with high H 2 S or CO 2 content. The phase transition near the wellhead could cause a great expansion of the gas volume, thus leading to catastrophic blowout if not controlled well, such as the CO 2 blowout at New Mexico, Colorado and Wyoming in the U.S. and the H 2 S blowout accident in well 16 H, Luojia in China. Therefore, the research on multiphase flow with sour gas invasion has become a key area in recent years.
According to the physical properties and phase transition of CO 2 , Skinner [4] proposed that Supercritical fluid could easily cause blowout. Yuan et al. [5] studied the influence of Supercritical phase behaviors on well control safety and quantitatively analyzed gas density and solubility distribution in the wellbore. Zhang et al. [6] analyzed the accident cause induced by phase transition in consideration of supercritical CO 2 and H 2 S phase behavior characteristics. Shi et al. [7] considered Z-factor, viscosity and density variations of H 2 S gas to establish a gas fluid flow model during H 2 S invasion, in which gas solubility in drilling fluid is ignored. Gao et al. [8] analyzed sour natural gas physical properties near the phase transition point by using a sample well in Puguang gas field, but the impact of wellbore heat transfer was not considered. Dou et al. [9] combined wellbore two-phase flow and steady-state heat transfer model based on phase transition and gas solubility to study the influence of sour gas density, pressure distribution and wellhead back pressure on phase transition position. However, he ignored cuttings and the physical properties variation of drilling fluid on HTHP. Sun et al. [10] established a multiphase flow model with considering the phase transition and the solubility of sour gas, while only H 2 S, instead of CO 2 component was analyzed, and the simulation of wellbore circulating temperature was not addressed. Yin et al. [11] proposed an improved approach to control the gas kick in high-pressure sour gas wells, however, it was not theoretically based on the multiphase flow model.
In deep HTHP formation, the physical parameters of the drilling fluids could change greatly along with the change of the depth. The wellbore temperature and pressure mutually influence each other, and the changes in the temperature and pressure will cause variations in the thermophysical and flow properties of the components of the wellbore fluids. However, most of the above studies ignored these factors, and none of them gave a quantitative comparison between H 2 S and CO 2 .
Thus, based on three governing equations (mass, momentum and energy conservation), a new multiphase flow model for sour gas kicks has been established with the consideration of sour gas solubility, phase transition and effects of temperature and pressure on the physical parameters of drilling fluid. And a transient heat transfer and pressure coupling solution is proposed to ensure the convergence and accuracy of the numerical simulation. On the basis of model validation, a case study is presented to illustrate the behavior of the wellbore temperature and pressure, and a parametric sensitivity analysis has been further conducted to provide a qualitative evaluation. These research results could provide some theoretical guidance for the exploration and development of deep sour gas reservoirs.

Fundamental Assumptions
To simplify the coupled wellbore flow and heat transfer mathematical model and its calculation, the following assumptions are given as follows: 1. Drilling fluid flows through the wellbore in one dimension, ignoring radial variation. 2. The drilling fluid is water-based, considering the dissolution and precipitation of sour gas in the drilling fluid. 3. The change of internal energy caused by phase transition of sour gas along the wellbore is not considered. 4. The heat transfer in the pipe wall, casing and formation are all in an unsteady state , thus, the axial heat conduction should be considered. 5. The axial heat conduction of drilling fluid in the wellbore could be ignored compared with the axial heat convection.

Mass Conservation Equation
Considering the infulence of the solubility of the sour gas, which could be constantly precipitated and dissolved along the wellbore with the change of temperature and pressure, the wellbore fluid components in this model are divided into four parts: drilling fluid, dissolved gas, free gas and cuttings.
where A is annular area, m 2 ; ρ l , ρ d , ρ f , ρ s are the density of the drilling fluid, the dissolved gas, the free gas and the cuttings respectively, kg/m 3 ; α l , α d , α f , α s are the volume fraction of the drilling fluid, the dissolved gas, the free gas, the cuttings respectively, dimensionless; v l , v d , v f , v s are the velocity of the drilling fluid, the dissolved gas, the free gas, the cuttings respectively, m/s; q g is gas invasion rate, kg/s; q d is gas dissolving rate, kg/s; q s is cuttings production rate, kg/s.

Pressure Drop Equation
According to Mass and Momentum Conservation Equations, the pressure drop of single-phase fluid flowing in drill pipe is given as follows.
where g is acceleration of gravity, 9.81 m/s 2 ; θ is deviation angle,°; f is fanning fraction factor, dimensionless; D d is drill pipe diameter, m.
Accurate flow pattern identification is an important prerequisite for establishing a comprehensive multiphase flow mathematical model [12]. According to previous research [13][14][15], two-phase flow patterns in wellbore are divided into five parts: bubble flow, dispersed bubble flow, slug flow, churn flow and annular flow. The pressure drop and drift flow coefficient along the well differ with the flow patterns. The specific algorithm and formula are not described here. A previous paper by He et al. [16] has introduced the algorithms and formulas in detail. Annular pressure drop of multiphase can be illustrated as follows.
where ρ m is gas-liquid-solid mixed fluid density, kg/m³; v m is gas-liquid-solid mixed fluid velocity, m/s; D h is annulus equivalent diameter, m.

Heat Transfer Model
During normal drilling fluid circulation, heat exchanges occur between formation and fluid in annulus and between fluid in drill pipe and annulus. The whole cycling process could be regarded as a heat exchanger with some boundary conditions. The physical model of wellbore-formation heat transfer is shown in Fig. 1. Based on the work by Abdelhafiz et al. [17], the heat transfer mathematical model is divided into five parts in radial direction: ① drill pipe, ② pipe wall, ③ annulus, ④ casing and ⑤ formation, and the corresponding formulas are as follows, respectively.
where ρ is density, kg/m 3 ; α is volume fraction, dimensionless; Q is the wellbore friction power, W/m; q is flow rate, L/s; T is temperature,°C; h is convective heat transfer coefficient, W/(m·°C); k is heat conductivity coefficient, W/(m 2 ·°C); C is specific heat, J/(kg·°C); z is the axial coordinate, m; r is the radial coordinate, m. The subscript: l is liquid phase; g is gas phase; p is drill pipe; w is pipe wall; a is annulus; c is casing; f is formation; pi is the inner wall of drill pipe; po is the outer wall of drill pipe; ci is the inner wall of casing; co is the outer wall of casing.
The convective heat transfer coefficient of gas-liquid fluid is quite different from that of single-phase flow. The convective heat transfer coefficient correlations proposed by Gao [18] is adopted in this paper: Aggour model for bubble flow, Knott model for dispersed bubble flow, Rezkallah and Sims model for slug and churn flow, Ravipudi and Gobold model for annular flow.

Sour Gas Z-Factor
There are many sour gas Z-factor calculating methods so far, which mainly can be divided into two categories: empirical correlation and state equation method. Common empirical correlations include DPR model, DAK model, HY model and HTP model. Among them, DPR and DAK model combining WA correction method are the most accurate, with less than 1.3% absolute average error [19]. However, these methods can be only applied within a certain range of temperature and pressure and may generate blind spots in calculation. Therefore, PR-EOS model [20], which is relatively accurate and suitable for the entire temperature and pressure range, is adopted to calculate sour gas Z-factor in this study Eq. (12) can be turned into a cubic equation of Z-factor: where p ci is the critical pressure of component i, Pa; T ci is the critical temperature of component i, K; R is universal gas constant, 0.008314 MPa·m 3 /(kmol·K); w i is the Z-factor of component i, dimensionless; T ri is reduced temperature of component i, dimensionless; k ij is the binary interaction coefficient [21], dimensionless; N is the number of gas composition.

Sour Gas Viscosity
Wu et al. [22] used empirical model, including Dempsey, Lee, LBC, DS and Londono model, as well as three correction formula of sour gas, YJS, Standing and Elsharkawy, to calculate the viscosity of sour gas under different temperature and pressure, with a total of 188 sets of data. The results showed that the average relative error could be minimum, only 3.4%, when using Dempsey model combining with Standing correction method.

Sour Gas Solubility
The solubility of H 2 S and CO 2 in water-based drilling fluid is much higher than that of CH 4 . As the gas moves upward along the wellbore, the solubility decreases with the decline of temperature and pressure, and the dissolved gas is gradually separated from the drilling fluid, which will have a great impact on the wellbore fluids properties. Duan et al. [23] established a solubility model of sour gas (CH 4 , H 2 S and CO 2 ) in aqueous solution by using the equation of state and the theory of interaction between specific particles. That model, widely used in the world, reproduced hundreds of sets of experimental data with high computational accuracy and an average relative error of no more than 7%. Therefore, this study adopted the model to calculate the solubility.

Drilling Fluid Density and Plastic Viscosity
Physical parameters of the drilling fluid, such as density and plastic viscosity vary with the change of the temperature and pressure. The following formulas show the relations between density and plastic viscosity of water-based drilling fluid and the temperature and pressure, which are determined by laboratory experiments.
where P 0 is surface pressure, MPa; T 0 is surface temperature, K; ρ(P,T) is drilling fluid density under pressure is P and temperature is T, kg/m 3 ; μ(P,T) is drilling fluid plastic viscosity under pressure is P and temperature is T, mPa·s; ρ 0 is drilling fluid density on the ground, kg/m 3 ; μ 0 is drilling fluid plastic viscosity on the ground, mPa·s.

Boundary Conditions
During the whole drilling process, the wellhead back pressure remains unchanged as the wellbore pressure boundary condition.
The drilling fluid temperature at the drill string inlet can be measured directly. The temperature of drilling fluid in the drill pipe, the pipe wall and the annulus are the same at the bottomhole.

Finite Difference Equations
A set of non-linear governing Eqs. (1)-(11) is constructed, including wellbore flow and heat transfer model. Considering the complicated partial differential form of the governing equations, we use a finite difference numerical method to solve it. The governing equations can be discretized by the following principles: i) The first order spatial/time derivative adopts the backward difference format; ii) The second order spatial derivative adopts the central difference format. Consequently, the finite difference equations of the heat transfer model are expressed as: Q p À Á n j À q l q l C l ð Þ n j T n 1;j À T n

1;jÀ1
Δz À 2pr pi h pi À Á n j T n 1;j À T n 2;j ¼ q l C l ð Þ n j pr 2 pi T n 1;j À T nÀ1 where PS1 ¼ q g C g q g a g þ q l C l q l a l þ q s C s q s a s ; PS2 ¼ q g C g a g þ q l C l a l þ q s C s a s ; i, j and n represent radial, axial and time nodes, respectively.

Solution Algorithm
The influence of the wellbore-formation transient heat transfer should be considered to ensure the convergency and accuracy of the model. The specific steps of the wellbore temperature and pressure coupling solving method are as follows: 1. Discrete time and space by meshing the wellbore. 2. Put the known wellbore temperature T i n at time n into pressure drop formulas in drill pipe and annulus to solve wellbore pressure p i n and thermophysical properties of each component.
3. Assume the temperature T i n+1(0) , pressure p i n+1(0) and heat physical properties of each component at time n + 1 according to the calculated value at time n, then calculate a new wellbore temperature T i n+1 at time n + 1 by using fully implicit Crank-Nicolson method [24]. 4. Calculate a new wellbore pressure p i n+1 by putting known parameters into pressure drop equation. If ||p i n+1 -p i n+1(0) ||<ε p and ||T i n+1 -T i n+1(0) ||<ε T , the predictions of p i n+1(0) and T i n+1(0) are accurate enough, otherwise return to Step (3) until the requirement of accuracy is met. 5. The calculation continues until pre-set time is reached, ending the wellbore pressure and temperature coupling calculation.

Model Validation
Large-diameter pipe flow experiments performed at a variety of phase flow rates and inclinations angle were introduced in detail by Shi et al. [25]. The flow loop was 10.9 m long, 0.152 m in diameter, and the pipe inclination varied from 0°(vertical) to 90°(horizontal). The tap water and nitrogen were chosen as the liquid and gas phase respectively. During the experiments, the water and gas entered the pipe and flowed along the test section, where the flow pattern could be observed visually. The average liquid volume fraction was measured by using shut-in method, and this parameter was employed for model validation, because it is one of the most important physical quantities to characterize multiphase flow in wellbores. Fig. 2 shows the contrast between the experimental liquid volume fraction and the simulated results. Three sets of data for different inclinations range from 0°to 80°are used, in which the inclination equals FDMP, 2020, vol.16, no.5 0°, 45°and 80°. It should be noted that dashed lines in the figure indicate errors of ±10% and ±20%. It is obvious that most data are within the range from 0 to 20%, which have a clear tendency toward underestimate. The absolute average error is 9.54%. Hence, the simulated results of the liquid volume fraction match well with the measured data.

Simulation Results Analysis
A vertical sample well with water-based drilling fluid was simulated in the case study. H 2 S and CO 2 (the rest is CH 4 ) of different contents were selected as the invading sour gas to do the contrastive analysis. Furthermore, a parametric sensitivity analysis has been conducted to evaluate qualitatively. The detailed input data of this simulation are shown in Tab. 1.   Fig. 3 shows the annulus temperature profiles with pure CH 4 , H 2 S and CO 2 invasion after 8 hours circulation. The results illustrate that the annulus temperature is significantly different from the formation temperature because of the wellbore-formation heat transfer. The highest temperature in circulating annulus exist at 1/10 of the depth above the bottomhole (nearly 4500 m) rather than at the bottomhole. However, the highest circulating temperature of single-phase drilling fluid in annulus usually exist at 1/6~1/7 of the depth above the bottom [26], which is mainly caused by the difference of heat transfer coefficient between gas-liquid two-phase and single-phase. According to the results shown in Figs. 4a and 4b, we observe that the outlet temperature increases gradually with the decrease of the bottomhole temperature, and both tend to be stable at last. Besides, during the circulation, the annulus temperature with H 2 S kicks is the highest, followed by CO 2 , and CH 4 the lowest.
As can be seen from Fig. 5, after 8 hours circulation, the annulus pressure with H 2 S kicks is higher than that with CH 4 kicks, and that with CO 2 invasion falls in between. The reason is that the solubility of H 2 S and CO 2 in the wellbore is much higher than that of CH 4 , most of which exist in the form of dissolution, and the free gas in the wellbore is greatly reduced. Moreover, the solubility of H 2 S is higher than that of CO 2 under the same conditions. The pressure differences among the three curves increase with the increase in the well depth, reaching maximum at the bottom. The bottomhole pressure difference between H 2 S and CO 2 invasion is 2.93 MPa, and that between CO 2 and CH 4 is 5.92 MPa. Fig. 6 shows the gas density profiles in the wellbore with CH 4 , H 2 S and CO 2 invasion after 8 hours circulation. The density of CH 4 varied within a small range in the wellbore compared with that of CO 2 and H 2 S. The H 2 S density along with the well has the highest rate of change. The density of CH 4 which is an ordinary hydrocarbon, varies relatively slow along the wellbore and no dramatical volume expansion exists. When considering the H 2 S or CO 2 sour gas invasion, the gas expansion is not evident in the early stage,but when the gas is moves to the phase transition point close to the wellhead, it turns to gaseous state from liquid or supercritical state, resulting in a dramatical volume expansion. This may cause instant blowout which almost leaves no time for the oilfield personnel to respond. The phase transition points of H 2 S and CO 2 are at 600 m and 1200 m depth respectively. The transition point of H 2 S is closer to the wellhead and the gas expands faster, which makes the gas kick more dangerous, imperceptible and the well control more difficult. Fig. 7 shows the free gas volume fraction profiles in the wellbore with CH 4 , H 2 S and CO 2 invasion after 8 hours circulation. With the same gas invasion rate, the CH 4 free gas volume fraction is the largest, CO 2 comes second and H 2 S is the least. Compared with CO 2 and CH 4 , the free gas volume fraction of H 2 S is with the highest change rate along with the change in well depth. The free phase volume fraction of H 2 S keeps zero below 2000 m depth, which means that almost all H 2 S is dissolved. The H 2 S and CO 2 are in supercritical state at the bottom, thus the volume expansion rate is very small compared with CH 4 , and little gas is released out. However, once the gas phase transition occurs due of the continuous drop of the temperature and pressure caused by the gas slipping up along the wellbore, the free gas volume fraction will dramatically rise. Thus, the sour gas invasion is difficult to detect at the early stage, and is also hard to control after detection.   8 shows the physical parameters profiles of the drilling fluid during circulation. t = 0 h refers to the original formation temperature, both the drilling fluid density and plastic viscosity gradually decrease with the increase of the well depth, indicating that the influence of temperature is dominant. Under the influence of the circulating temperature field, both the drilling fluid density and plastic viscosity first decrease and then increase with the increase of the well depth. This is the result of the dual action of temperature and pressure on the drilling fluid. The physical parameters (density, plastic viscosity) of drilling fluid are greatly affected by the circulating temperature. The variation in the density and plastic viscosity of drilling fluid needs to be fully considered when establishing the hydrodynamics model of HTHP, so as to accurately calculate the wellbore pressure and temperature and to provide reliable parameters for drilling operation.

Sour Gas Invasion with Different H 2 S/CO 2 Content
Figs. 9 and 10 show the sour gas solubility profiles with different H 2 S/CO 2 content after 8 hours circulation. The H 2 S and CO 2 in supercritical state are characterized by low viscosity, large diffusivity and high solubility. Under the bottomhole condition, a large amount of sour gas dissolves into the waterbased drilling fluid in a supercritical state. In the gas slipping up process, the solubility decreases and the amount of dissolved gas reduces gradually with the drop of the temperature and pressure. The amount of dissolved gas decreases slowly before reaching the phase transition point, afterward it decreases significantly. When the content of H 2 S/CO 2 increases, the amount of dissolved gas increases, and so does its decreasing rate. Under the same condition, the sour gas with high H 2 S content is with larger solubility and higher change rate than that with CO 2 content.  The simulation results of bottomhole pressure with different H 2 S/CO 2 content after 8 hours circulation are shown in Fig. 11. The bottomhole pressure increases linearly with the content of H 2 S/CO 2 increases. With the same amount of sour gas, the bottomhole pressure with H 2 S kicks is higher than that with CO 2 . The pressure difference gradually increases with the increase of H 2 S/CO 2 content. When the sour gas content is 10%, the pressure difference is only 0.12 MPa. However, when the sour gas content reaches 70%, the pressure difference comes up to 1.75 MPa.

Parametric Sensitivity Analysis
Bottomhole temperature and pressure are among the most concerned issues in the drilling operation, which are often influenced by many variables. Therefore, a parametric sensitivity analysis has been conducted to quantitatively evaluate the effects of the selected nine parameters on the bottomhole temperature and pressure. The parameter values in Tab. 1 are used for inputs, pure CH 4 is simulated and the circulation time is 8 hours. The results of sensitivity analysis are shown in Tab. 2. It should be noted that the conclusions drawn from the sensitivity analysis can be only applied under the specific conditions used in the simulation.
The geothermal gradient is an important parameter affecting the temperature. Considering the difference between the values of the lower and the upper geothermal gradient reaches 33.3%, the bottomhole temperature changed 26.1%, and the pressure changed 2.5%.
Bottomhole temperature and pressure increase linearly with the increase of the depth. The difference (40%) in well depth results in an ultimate change of 52.8% and 41.1% in the bottomhole temperature and pressure respectively. The bottomhole pressure increases linearly with the increase of the back pressure, but the effect of the back pressure on the bottomhole temperature is negligible.
The difference (40%) between the lower and the upper gas kick rate results in an ultimate change of 3.7% and 7.3% in the bottomhole temperature and pressure respectively. While the difference (50%) between the lower and the upper pump rate results in an ultimate change of 5.6% and 15.5% in the bottomhole temperature and pressure respectively.  The drilling fluid density ρ 0 and plastic viscosity u 0 under surface conditions are selected as independent variables to analyze the effect of the physical parameters of the drilling fluid on the bottomhole temperature and pressure. The difference between the lower and the upper drilling fluid density (16.7%) caused an ultimate change of 3.5% and 18.6% in the bottomhole temperature and pressure respectively, while the effects of the plastic viscosity on the bottomhole temperature and pressure are negligible.
In conclusion, the results of the parametric sensitivity analysis are given as follows: i) As for the effect on bottomhole temperature, depth > geothermal gradient > drilling fluid density > pump rate > gas kick rate > back pressure > drilling fluid plastic viscosity; ii) As for the effect on bottomhole pressure, drilling fluid density > depth > pump rate > gas kick rate > geothermal gradient > back pressure > drilling fluid plastic viscosity; iii) The effects of H 2 S content on bottomhole temperature and pressure are greater than that of CO 2 content.

Conclusions
1. Considering the sour gas solubility, phase transition and the effects of temperature and pressure on the physical parameters of drilling fluid, a new multiphase flow model with transient heat transfer and pressure coupling for sour gas invasion has been established and proved to be more accurate compared with the large-diameter pipe flow experimental data. 2. The position with the highest temperature in cycling annulus is at about 1/10 of the depth above the well bottom. During the circulation, the annulus temperature with H 2 S kicks is the highest, followed by CO 2 , and CH 4 is the lowest. It is difficult to detect the sour gas invasion at the early stage, and is also hard to control after detection. The phase transition point of H 2 S is closer to wellhead compared with CO 2 , resulting in a faster expansion rate, which is more imperceptible and dangerous. 3. Under the influence of the circulating temperature field, the drilling fluid density and plastic viscosity both first decrease and then increase with the increase of the well depth. With the same amount of sour gas, the bottomhole pressure is higher when H 2 S kicks than CO 2 kicks, and the pressure difference gradually increases with the increase of H 2 S/CO 2 content. 4. For the effect on bottomhole temperature, depth > geothermal gradient > drilling fluid density > pump rate > gas kick rate > back pressure > drilling fluid plastic viscosity. While as for the effect on bottomhole pressure, drilling fluid density > depth > pump rate > gas kick rate > geothermal gradient > back pressure > drilling fluid plastic viscosity. The effects of H 2 S content on bottomhole temperature and pressure are greater than those of CO 2 content.