A New Model of Bubble Migration Velocity in Deep Water Wellbore Considering Hydrate Phase Transition

1. Introduction

There are abundant oil and gas resources in deepwater area of the world, and the scale of oil and gas exploration and development is gradually increasing. However, the harsh marine environment and complex geological condition greatly increase the difficulty of oil and gas exploration and development [1]. The well will be shut-in when encountering severe sea conditions or complex conditions such as gas invasion during drilling. Gas in the reservoir invades the wellbore annulus through diffusion and displacement during well shut-in. When the diameter of the bubble invading the annulus is larger than the critical suspension diameter, the bubble will slide off and rise [2]. As the fluid temperature in the wellbore annulus is close to the environment temperature during well shut-in, when the bubble slips up to the hydrate generation area, the hydrate shell will be formed on the surface and continuously thicken. When the gas bubble coating the hydrate shell (hereinafter referred to as "hydrated bubble") migrate to the wellhead of the submarine mudline, the hydrate shell continues to grow and hydrated bubble accumulate under the condition of low temperature and high pressure. This can even cause serious accidents such as blowout preventer blockage, which will bring severe challenges to subsequent well opening operations [3]. Therefore, accurate prediction of bubble migration velocity is crucial for calculating the safe shut-in period of deepwater oil and gas wells.

The migration behavior of bubble in deepwater wellbore is closely related to the bubble shape, fluid properties and other factors. Due to the special temperature and pressure environment, when bubble enter the hydrate formation area, the influence of hydrate phase transition on the migration behavior of bubble is more complex [4]. In recent years, many scholars have studied the migration behavior of bubble in deepwater wellbore. Joseph et al. [5] established a model for the rising velocity of uncovered hydrate shell bubble (hereinafter referred to as "clean bubble") in deepwater environment, considering the influence of surface tension and fluid viscosity. Peebles and Garberd et al. [6] considered the influence of bubble deformation and flow field on the migration velocity of bubble, and established a segmented prediction model for the migration velocity of clean bubble based on the Reynolds number range. Collins et al. [7] further considered the influence of wall limitation on the rising velocity of deepwater bubbles, and modified the segmented prediction model established by Peebles and Garberd. Margaritis et al. [8] obtained a correlation formula for the migration velocity of clean bubble using a large number of experimental data on bubble migration, which is suitable for various non-Newtonian fluids. Taylor et al. [9] used a microscope to accurately measure the thickness change of the hydrate shell at the gas-liquid interface during the experiment. Bigalke et al. [10] studied the influence of shape parameters such as bubble aspect ratio on the migration velocity of methane and carbon dioxide bubbles during their ascent in deepwater environment. Sun et al. [11] established a dynamic model for the inward growth of hydrate shell on the surface of bubble in deepwater wellbore considering the permeation and mass transfer of water. Liu et al. [12] established a thickness prediction model for bubble migration in deepwater wellbore by considering factors such as internal pore renewal and mass transfer during hydrate shell growth. Sun et al. [13] explored the influence of undercooling on the thickness of hydrate shell, and obtained the undercooling influence coefficient by fitting the experimental data. Wei et al. [14] established a model for the migration velocity of hydrate bubble based on the force balance in the hydrate formation area, and characterized the increase of bubble density caused by the hydrate shell coating behavior with equivalent density.

The investigation shows that previous studies have been carried out on the prediction of the migration velocity of clean bubble, the experimental observation and model prediction of the growth thickness of the hydrate shell. Since the migration of hydrated bubble involves complex processes such as gas-water bidirectional mass transfer, pore renewal in hydrate shell, gas diffusion and dissolution, and gas properties changing with temperature and pressure, there are few studies on the migration velocity of hydrated bubble in deepwater wellbore environment. In addition, the prediction model of bubble migration velocity in the whole process from generation to migration to the subsea wellhead is rarely studied. However, the research on the velocity of bubble migration in the whole process from generation to migration to the subsea wellhead is very important for accurately predicting the safe shut-in period and subsequent safe well-opening. Therefore, based on the characteristics of bubble migration in deepwater wellbore, this study established a velocity prediction model for the whole process of bubble migration from bottom hole to underwater wellhead by considering factors such as mass transfer, phase transition, temperature and pressure change, diffusion and dissolution. The model can accurately characterize the complex processes such as the dynamic change of hydrate shell thickness, gas-liquid bidirectional mass transfer and pore renewal, and the accuracy of the model is verified by experimental observation data. Using the data of a deepwater well in the South China Sea, the influence of the initial size of the bubble and the properties of the annular fluid on the migration velocity were analyzed. It is hoped that this study can provide a reference for quantitative description and characterization of complex bubble migration behavior with phase transition and mass transfer in deepwater environment.

2. Model description

Reservoir gas invades the wellbore annulus through diffusion and displacement from the rock pores in the open-hole section during the shut-in period. The gas invading the wellbore is dispersed in the annulus fluid in the form of bubbles. When the bubble size is larger than the critical suspension condition, it begins to slip and rise. Sikorskir et al. found in their experiments that bubble width can better reflect the influence of fluid resistance on bubble suspension than equivalent diameter, and proposed the critical suspension condition considering the maximum bubble width [15]:

$$Y_g=\frac{2\pi {\tau }_y{R}_{\max }^2}{\Delta \rho g{V}_b}$$

(1)

where, $Y_g$ is the dimensionless yield stress. $\Delta \rho$ is the gas-liquid two-phase density difference, kg/m³. ${\tau }_y$ is yield stress, Pa. $g$ is the acceleration of gravity, m/s². $V_b$ is bubble volume, m³. $R_{\max }$ is the maximum width of the bubble in the horizontal direction, m.

When the maximum bubble width is larger than the critical suspension width, it migrates upward in the annulus fluid, and its migration process is shown in Figure 1. In Figure 1, the blue curve is the phase equilibrium curve of methane hydrate, and the orange curve is the shut-in wellbore temperature curve. The corresponding depth range of the area surrounded by the two curves represents the hydrate formation area. When bubble migrates to this area, they change from clean bubble to hydrated bubble. When the bubble is below the hydrate formation area, its migration process is mainly affected by temperature, pressure and the properties of the drilling fluid. The migration process of gas bubble in the hydrate formation area is also affected by the degree of supercooling and mass transfer rate. The final position of bubble migration during well shut-in is the wellhead of submarine mudline. Therefore, this study mainly establishes the velocity model of the whole process of bubble migration from bottom hole to underwater wellhead, including the clean bubble section and the hydrated bubble section.

Bubble migration in wellbore annulus is accompanied by the comprehensive effect of various forces, as shown in Figure 2. These include the gravity and buoyancy force of the bubble, the viscous resistance caused by the apparent viscosity of the fluid to the rise of the bubble, the additional mass force caused by the acceleration of the bubble and the movement of the surrounding fluid, and the Bassett force caused by the instantaneous velocity change of the bubble. Compared with the clean bubble, the number of forces on the hydrated bubble has not changed, but the force magnitude and change rate are different. Therefore, the clean bubble section and the hydrated bubble section can be unified modeling.

The gravity and buoyancy of bubble can be expressed as:

$$F_g=\frac{4}{3}\pi g\left[{\left(r_g+\delta \right)}^3{\rho }_h-{r_g}^3\left({\rho }_h-{\rho }_g\right)\right]$$

(2)

$$F_v=\frac{4}{3}\pi {\left(r_g+\delta \right)}^3{\rho }_{l}g$$

(3)

where, $F_g$ is the force of gravity on the bubble, N. $r_g$ is bubble equivalent radius, m. ${\rho }_g$ and ${\rho }_l$ are the density of bubble and liquid, kg/m³. ${\rho }_h$ is the density of hydrate shell, kg/m³. $\delta$ represents the thickness of hydrate shell, m. $F_v$ is the buoyancy force of the bubble, N.

The viscous resistance caused by bubble migration in non-Newtonian fluid is mainly controlled by liquid viscosity, bubble volume and bubble migration velocity. The viscous resistance of the bubble can be expressed as [16]:

$$F_D=\frac{1}{2}{C}_D\pi {\rho }_l{\left(r_g+\delta \right)}^2{v}_g^2$$

(4)

where, $F_D$ is the viscous resistance of the bubble, N. $C_D$ is the drag coefficient, dimensionless. $v_g$ is the migration velocity of bubbles, m/s.

As the bubble density is far less than the fluid density in the wellbore annulus, the bubble will accelerate when it rises from the bottom of the well, and drive the fluid around the bubble to move. The resistance increased by driving the fluid motion is the additional mass force, and its expression is [17]:

$$F_m=\frac{2}{3}\pi {\left(r_g+\delta \right)}^3{\rho }_l\frac{d{v}_g}{dt}$$

(5)

where, $F_m$ is the additional mass force on the bubble, N.

Bubbles will deform during the rising process, which will generate instantaneous flow resistance, whose expression is [18]:

$$F_b=K_b{\left(r_g+\delta \right)}^2\sqrt{\pi \mu {\rho }_l}{\displaystyle {\int }_{t_0}^t\frac{d{v}_g/d{t}_{\lambda }}{\sqrt{t-t_{\lambda }}}}d{t}_{\lambda }$$

(6)

where, $F_b$ is Bassett force, N. $K_b$ is Bassett force coefficient. $\mu$ is fluid viscosity, Pa·s. $t_0$ is the initial time of deformation, s. $t$ is the deformation termination time, s. $t_{\lambda }$ is the integral variable, s.

Considering the effects of various forces in the process of bubble migration, the force relationship is as follow:

$$m_g\frac{d{v}_g}{dt}=F_v-F_g-F_d-F_m-F_b$$

(7)

where, $m_g$ is the mass of bubble, kg.

Substituting Equations (2) to (6) into Equation (7), the governing equation of bubble migration can be obtained as follow:

(8)

Equation (8) is the governing equation for the migration velocity of hydrated bubble in the deepwater wellbore annulus. When the bubble is in the non-hydrate formation region, the governing equation for the migration velocity of clean bubble can be obtained by substituting $\delta =0,{\rho }_h={\rho }_g$ into Equation (8):

$$\frac{4}{3}\pi r_g{\rho }_g\frac{d{v}_g}{dt}=\frac{4}{3}\pi g{r}_g({\rho }_l-{\rho }_g)-\frac{1}{2}{C}_D\pi {\rho }_l{v}_g^2-\frac{2}{3}\pi r_g{\rho }_l\frac{d{v}_g}{dt}-K_b\sqrt{\pi \mu {\rho }_l}{\displaystyle {\int }_{t_0}^t\frac{d{v}_g-d{t}_{\lambda }}{\sqrt{t-t_{\lambda }}}d{t}_{\lambda }}$$

(9)

3. Determination of key model parameters

According to Equations (8) and (9), the key parameters determining the migration velocity of hydrated bubble are bubble equivalent radius $r_g$, hydrate shell thickness $\delta$ and drag coefficient $C_D$, while the key parameters determining the migration velocity of clean bubble are equivalent radius $r_g$ and drag coefficient $C_D$. Since these key parameters are affected by temperature and pressure changes, mass transfer, and hydrate phase transition during the bubble migration process, the precise characterization of the key parameters in different bubble migration stages is the basis for accurately calculating the bubble migration velocity.

3.1. Dynamic growth thickness of hydrate shell

When bubble enters the hydrate formation area, a thin hydrate shell will rapidly form on its surface, and the thickness of the hydrate shell is usually called the initial thickness. Shi et al. [19] found in their experiments that hydrate nucleates on the surface of bubble, and then rapidly covers the whole bubble through lateral growth. Experimental statistics show that the formation time of the initial hydrate shell is between 10 and 25 s. Determining the initial hydrate shell thickness is the basis of establishing the hydrate shell thickness prediction model. Li et al. [20] measured the initial thickness of hydrate shell under different undercooling degrees. The results show that when the undercooling degree is greater than 1.0 K, the initial thickness of hydrate shell is inversely proportional to the undercooling degree.

$${\delta }_0=\frac{b}{\Delta T}$$

(10)

where, ${\delta }_0$ is the initial thickness of hydrate shell, m. $b$ is regression coefficient, m·K. $\Delta T$ is the degree of undercooling, K. According to the experimental data of Li et al., $b=3.1205\times {10}^{-5}$.

Lee et al. [21] observed the hydrate shell growing on the surface of water droplets in a gas-dominated system, and found that the hydrate shell showed inward growth phenomenon. Li et al. [22] experimentally studied the growth process of hydrate film on the surface of water droplets suspended in the oil phase, and also found the inward growth phenomenon of hydrate shell. Liu et al. [12] observed the inward growth phenomenon and local slight bulge of hydrate shell on the surface of suspended bubble in the aqueous phase dominated system. These studies indicate that gas-water mass transfer determines the growth direction and thickness variation of hydrate shell. During the dynamic growth of hydrate shell, the changes of temperature and pressure will lead to the formation of micro-cracks in its interior [22]. Water penetrates into the hydrate shell under the action of capillary force, while gas diffuses outward through the hydrate shell. As shown in Figure 3, water penetrates into the gas-side hydrate formation layer through the micro-cracks in the hydrate shell. The water molecules infiltrated into the hydrate shell will combine with the gas molecules diffused outward to form hydrate, which makes the internal structure of the hydrate shell compact. When the gas molecules enter the water-side hydrate formation layer by diffusion, they will first dissolve into the liquid phase. When the concentration of liquid gas reaches saturation, the excess gas molecules will accelerate the growth of hydrate on the water side, and when the concentration of liquid gas is lower than the saturation concentration, the dissolution of gas in the liquid phase will inhibit the growth of hydrate on the water side. Based on the above analysis, the dynamic growth of hydrate shell mainly includes outward growth, inward growth and internal pore renewal.

Based on the above analysis, the reasonable assumptions before establishing the dynamic model of hydrate shell thickness in this study are as follows:

(1) Hydrate shell is a kind of plastic-like material, which has a certain ability to resist damage [23]. Therefore, the mechanical failure of hydrate shell during bubble migration is not considered.

(2) The hydrate shell has a porous medium-like structure, consisting of hydrate crystals and microscopic pore throats.

After the formation of the initial hydrate shell, the subsequent thickness change depends on the hydrate formation rate inside and outside the hydrate shell. The thickness change of hydrate shell can be expressed as:

$$\frac{4}{3}\pi \frac{d\left\{\left[{\left(r_g+\delta \right)}^3-r_g^3\right]{\rho }_h\right\}}{dt}=h_{fi}-h_{fo}$$

(11)

where, $t$ is time, s. $h_{fi}$ and $h_{fo}$ are hydrate formation rates on the inside and outside of the hydrate shell respectively, kg/s.

The bubble equivalent radius $r_g$, hydrate shell thickness $\delta$ and hydrate shell density ${\rho }_h$ all change with time due to the influence of temperature and pressure change, mass transfer and hydrate shell compact growth. Therefore, the governing equation of hydrate shell thickness can be obtained by differentiating all variables in Equation (11):

$$\frac{d\delta }{dt}=\frac{1}{4\pi {\rho }_h{\left(r_g+\delta \right)}^2}\left(h_{fi}-h_{fo}\right)-\frac{2\delta }{r_g+\delta }\frac{d{r}_g}{dt}-\frac{\left[{\left(r_g+\delta \right)}^3-{r_g}^3\right]}{3{\rho }_h{\left(r_g+\delta \right)}^2}\frac{d{\rho }_h}{dt}$$

(12)

When the water molecules on the outside of the hydrate shell permeate to the inside of the hydrate shell through micro-cracks, a large number of gas molecules around it rapidly combine with water to form the hydrate. Therefore, the water penetration rate determines the thickening rate of the inner hydrate shell, and the inward hydrate formation rate can be expressed as:

$$h_{fi}=\frac{M_g+n{M}_w}{n{M}_w}{\rho }_w{f}_w$$

(13)

where, $n$ is the hydration number, dimensionless. $M_g$ is the molar mass of gas, kg/mol. $M_w$ is the molar mass of water, kg/mol. ${\rho }_w$ is the density of water, kg/m³. $f_w$ is the permeation rate of water through hydrate shell, m³/s.

According to the model proposed by Mori et al., the total mass transfer rate of water penetration can be expressed as follows [24]:

$$f_w=\frac{\pi }{8{\mu }_l\delta }\left(2\sigma +\frac{r_c\Delta p}{\cos \beta }\right){\gamma }_c$$

(14)

where, $\sigma$ is the gas-water interfacial tension, N/m. $\beta$ is the contact angle of capillary on the water side, rad. $r_c$ is the capillary radius, m. $\Delta p$ is the pressure difference between the inside and outside of the hydrate shell, Pa. ${\gamma }_c$ is the ability of water molecules to pass through micro-cracks, m³, which depends on the number and size of micro-cracks in the hydrate shell:

$${\gamma }_c=4\pi r_g^2{n}_d\frac{r_c^3\cos \beta }{s}$$

(15)

where, $n_d$ is the number of micro-cracks per unit area in the hydrate shell, 1/m². $s$ is the tortuosity of the micro-cracks, dimensionless.

On the outside of the hydrate shell, the formation rate of hydrate is mainly controlled by the mass transfer and dissolution rate of gas molecules. Therefore, the outer formation rate of hydrate shell can be expressed as:

$$h_{fo}=\frac{M_g+n{M}_w}{M_g}{\rho }_g\left(f_s-f_d\right)$$

(16)

where, $f_s$ is the mass transfer rate of inner gas diffusion through the hydrate shell, m³/s. $f_d$ is the dissolution rate of gas in the liquid phase, m³/s.

The mass transfer rate of bubble diffusion can be obtained by solving the steady-state diffusion equation. The diffusion equation adopted in this study is [25]:

$$d_h\left(\frac{{\partial }^2{C}_g}{\partial r}+\frac{2}{r}\frac{\partial C_g}{\partial r}\right)=0$$

(17)

where, $C_g$ is the gas concentration, m³/ m³. $d_h$ is the gas diffusion coefficient inside the hydrate shell, m²/s. The diffusion coefficient $d_h$ inside the hydrate shell reflects the diffusivity of the gas along the concentration gradient. Ogasawara et al. [26] obtained that the gas diffusion coefficient inside the hydrate is about 10^-11~10^-12 m²/s on the basis of experimental research.

In this study, the boundary conditions corresponding to the diffusion equation are:

$$r=r_g,C_g=C_{i}r=r_g+\delta ,C_g=C_o$$

(18)

where, $C_i$ is the gas concentration inside the hydrate shell, m³/m³. $C_o$ is the gas concentration outside the hydrate shell, m³/m³.

According to the boundary conditions of Equation (18), the analytical formula of Equation (17) is obtained as follow:

$$C_g\left(r\right)=\left[r_g{C}_i-r_g\frac{\left(r_g+\delta \right)C_o-r_g{C}_i}{\delta }\right]\frac{1}{r}+\frac{\left(r_g+\delta \right)C_o-r_g{C}_i}{\delta }$$

(19)

According to equation (19), the gas concentration gradient at position $r=r_g$ is:

$$\frac{\partial C_g}{\partial r}{|}_{r=r_g}=\frac{\left(r_g+\delta \right)}{r_g\delta }\left(C_o-C_i\right)$$

(20)

The mass transfer rate of inner gas diffusion through the hydrate shell is:

$$f_d=-4\pi r_g^2{d}_h\frac{\partial C_g}{\partial r}{|}_{r=r_g}=4\pi d_h\frac{r_g\left(r_g+\delta \right)}{\delta }\left(C_i-C_o\right)$$

(21)

The dissolution rate of gas in liquid phase is closely related to the gas concentration difference and mass transfer coefficient at the interface between hydrate shell and liquid phase. In addition, the rate of gas dissolution in the liquid phase also depends on the bubble surface area. The gas dissolution rate at the hydrate shell-liquid interface can be expressed as follow [27]:

$$f_r=4\pi k_{hl}{\left(r_g+\delta \right)}^2\left(C_h-C_l\right)$$

(22)

where, $C_h$ is the gas concentration at the hydrate shell-liquid interface, m³/m³. $C_l$ is the gas concentration in the liquid phase, m³/m³. $k_{hl}$ is the mass transfer coefficient between hydrate shell and liquid phase, m/s. The mass transfer coefficient determines the dissolution rate of gas in liquid phase. Rehder et al. [28] showed that the existence of hydrate shells would increase the mass transfer resistance of gas dissolution and reduce the mass transfer coefficient. The formation of the hydrate shell leads to the transition from the original mass transfer across the gas-liquid interface to the mass transfer across the solid hydrate shell, and the mass transfer coefficient of gas in the solid phase is significantly lower than that in the liquid and gas phases. In order to accurately characterize the effect of the dissolution rate of gas mass transfer across the hydrate shell on the thickness of the hydrate shell, error analysis and optimization of the typical mass transfer model of bubble dissolution in Table 1 should be carried out by using experimental data. Rehder et al. [28] experimentally studied the variation law of the equivalent radius of methane bubbles released in the water depth range from 400 m to 1500 m with the bubble migration time. The experimental data were collected and recorded by ROV (Remote Operated Vehicle) observation. The data are of high practical value and have a strong reference significance for the establishment and verification of bubble migration model. Therefore, in this study, error analysis and optimization of the typical mass transfer model of bubble dissolution were carried out based on the experimental data of bubble equivalent radius variation with migration time in hydrate shell formation segment obtained by Rehder et al., which released methane bubbles at a depth of 1098.0 m.

Figure 4 shows the comparison between the calculated results of different typical mass transfer models of bubble dissolution and the experimental data. In this experiment, the equivalent radius of the hydration bubble was recorded by ROV several times from 180 s to 470 s after bubble release. The Clift model, Levich model and Acrivos model all underestimate the gas mass transfer coefficient across the hydrate shell. Compared with the Rehder et al. experimental data, the maximum errors of the models are 21.26%, 17.35% and 23.61% respectively. The Oellrich model, Leclair model, Johnson model and Winnikow model all overestimate the gas mass transfer coefficient across the hydrate shell. Compared with the Rehder et al. experimental data, the maximum error of these models are 22.09%, 68.53%, 28.65% and 43.11% respectively. Based on the above analysis, the existing models can not accurately describe the gas dissolution mass transfer behavior of hydrated bubble. Among them, the average error of the Levich model is the smallest, because the situation that the gas-liquid interface is an immovable interface is considered in the modeling process. The reason for the error of this model is that the influence of the flow field around the bubble on the mass transfer coefficient is not considered. Therefore, a new mass transfer coefficient correlation model is obtained by considering the introduction of Reynolds number to modify the Levich model. The coefficients in the new model are regressed according to the experimental data obtained by Rehder et al.

$$k_{hl}=0.624{v_g}^{1/3}{\left(\frac{E}{r_g+\delta }\right)}^{2/3}\left({\mathrm{Re}}^{1/4}+0.72\right)$$

(23)

Error analysis of the new mass transfer coefficient model was carried out using experimental data obtained by Rehder et al. [28] releasing methane bubbles at water depth of 1209.6 m and 1511.4 m. In these two experiments, the equivalent radius of the hydrated bubbles was observed and recorded respectively during the 140 s to 900 s and the 80 s to 1100 s after the bubbles were released. Figure 5(a) and Figure 5(b) show the comparison between the calculated results of the new mass transfer coefficient model for bubble dissolution and the experimental data. The average error is 5.41% and 6.92% respectively, which indicates that the new mass transfer coefficient model of bubble dissolution can accurately characterize the mass transfer behavior of gas across hydrate shell.

The bidirectional mass transfer of gas and water across the hydrate shell is also accompanied by the renewal of pores in the hydrate shell. According to the relationship between the gas and water consumption during the formation of hydrate, we can obtain:

$$4\pi r_g^2{n}_d\left(\pi r_{c1}^{2}s\delta -\pi r_c^{2}s\delta \right)\frac{{\rho }_h}{M_g+n{M}_w}=n_{gc}$$

(24)

where, $r_c$ and $r_{c1}$ are pore radius before and after pore renewal respectively, m. $n_{gc}$ is the gas consumption, mol. The variation of pore radius in hydrate shell with time can be obtained from equation (24):

$$\frac{d{r}_c}{dt}=-\frac{M_g+n{M}_w}{8{\pi }^2{\rho }_{h}s{n}_d{r}_c{r}_g^2\delta }{m}_{gs}$$

(25)

where, $m_{gs}$ is the gas consumption rate in the process of pore renewal, mol/s. The difference between the gas concentration in the hydrate shell near the gas phase and the gas concentration in the hydrate shell near the liquid phase represents the speed of this consumption rate, which can be expressed as:

$$m_{gs}=4\pi r_g^2{d}_h\left(\partial C_g/\partial r\right){|}_{r=r_g}-4\pi {(r_g+\delta )}^2{d}_h\left(\partial C_g/\partial r\right){|}_{r=r_g+\delta }$$

(26)

3.2. Bubble dynamic equivalent radius

When the bubble is in the non-hydrate formation area in the wellbore annulus, the temperature and pressure change and dissolution are the main control factors affecting the equivalent radius of bubble. Dissolution refers to the effect of the bubble spreading into the surrounding liquid phase due to the difference in concentration during the migration process, which causes the bubble to shrink. In the process of bubble rising, the variation of temperature and pressure in deepwater wellbore annulus will lead to bubble expansion. Considering the influence of the above factors, the governing equation of equivalent radius in clean bubble migration section is established:

$$4\pi r_g^2{\rho }_g\frac{d{r}_g}{dt}=\frac{4}{3}\pi r_g^3{v}_g\left({\rho }_{l}g\frac{\partial {\rho }_g}{\partial P_a}+T_g\frac{\partial {\rho }_g}{\partial T_a}\right)-4\pi r_g^2{k}_{gl}\left(m_{g-l}-m_l\right)$$

(27)

where, $P_a$ is the annular pressure of deepwater wellbore, Pa. $T_a$ is deepwater wellbore annulus temperature, K. $T_g$ is the temperature gradient of annulus drilling fluid, K/m. $k_{gl}$ is the mass transfer coefficient at the gas-liquid interface, m/s. $m_{g-l}$ is the gas mass concentration at the gas-liquid interface, kg/m³. $m_l$ is the gas mass concentration in drilling fluid, kg/m³.

In Equation (27), the first term on the right represents the influence of wellbore annular temperature and pressure variation on the bubble equivalent radius, and the second term on the right represents the influence of dissolution on the bubble equivalent radius. The mass transfer coefficient determines the effect of dissolution on bubble shrinkage. Clift et al. conducted relevant research on the dissolution and mass transfer of clean bubble, and established multiple mass transfer coefficient models for different bubble migration states. Considering the air bubble migration process in deepwater wellbore annulus, this study adopts the clean bubble dissolution mass transfer coefficient model at low Reynolds number in deepwater established by Clift et al. [29]:

$$k_{gl}=\frac{2}{\sqrt{\pi }}\left[1-\frac{2}{3}\frac{1}{{(1+0.1415{\mathrm{Re}}^{2/3})}^{3/4}}\right]{\left(\frac{E{v}_g}{2r_g}\right)}^{1/2}$$

(28)

where, $\mathrm{Re}$ is the Reynolds number, dimensionless. $E$ is the diffusion coefficient of gas in drilling fluid, m²/s.

When bubble migrate to the hydrate formation area, the change of bubble equivalent radius caused by the inward growth of hydrate shell can be described by the following equation:

$$\frac{M_g+n{M}_w}{n{M}_w}{\rho }_w{f}_w=-4\pi r_g^2{\rho }_h\frac{d{r}_g}{dt}$$

(29)

3.3. Drag coefficient during bubble migration

In the research on resistance coefficient of bubble migration, the existing models are mainly established by correlating with dimensionless numbers such as Reynolds number and Morton number, and the coefficients in the correlation model are obtained by fitting a large number of experimental data. The drag coefficient model of the clean bubble migration section is relatively mature. Since this paper mainly studies the migration of bubble in non-Newtonian fluid in deepwater wellbore annulus, the drag coefficient model established by Rodrigue et al. [36] based on the experimental data of bubble migration in non-Newtonian fluid is adopted.

$$C_D=\frac{16}{\mathrm{Re}}\left[2^{n_w-1}{3}^{\frac{n_w-1}{2}}\frac{1+7n_w-5{n_w}^2}{n_w\left(n_w+2\right)}\right]$$

(30)

where, $n_w$ is the fluidity index, dimensionless.

Since the current research on the drag coefficient in the migration section of hydrated bubble is not mature, this study selects the drag coefficient model suitable for the migration of hydrated bubble by comparing the previous models with experimental data. As the applicable conditions of each model are different, the drag coefficient model which is widely used to describe bubble migration is selected for comparative analysis with the experimental data obtained during the previous research of our team [37]. By comparing the agreement between the calculated results of each model and the experimental data, the drag coefficient model suitable for the hydrated bubble migration section in non-Newtonian fluid is optimized. The typical drag coefficient model to describe bubble migration is shown in Table 2.

As shown in Figure 6(a) and Figure 6(b), the Bigalke model has a small difference with the experimental value at $0.5<C_D<1.2$, the Tomiyama model has a high coincidence with the experimental value at $C_D>1.5$, and the Bozzano model has a high coincidence with the experimental value at $0.95<C_D<1.75$. The average error of these three models in the above interval is less than 10.0%. Since the drag coefficient calculation model is usually selected according to the Reynolds number of bubble migration, the relationship between Reynolds number and drag coefficient in the process of hydrated bubble migration obtained from previous experiments [37] is used in this study to divide the applicable range of Bigalke model, Tomiyama model and Bozzano model.

As shown in Figure 7, when $\mathrm{Re}<275$, the drag coefficient value is roughly between 0.5 and 0.95, and the Bigalke model has good adaptability. When $275\le \mathrm{Re}\le 445$, the drag coefficient is roughly between 0.95 and 1.5, and the Bozzano model has good adaptability. When $\mathrm{Re}>445$, the drag coefficient is greater than 1.5, and the Tomiyama model has good adaptability.

$$C_D=\left\{\begin{array}{l}f{\left(\frac{a}{R}\right)}^2,\mathrm{Re}<275\\ F{\left(\frac{A}{Ro}\right)}^2,275\le \mathrm{Re}\le 445\\ \max \left\{\frac{24}{\mathrm{Re}}\left(1+0.15{\mathrm{Re}}^{0.687}\right),\frac{8}{3}\frac{Eo}{Eo+4}\right\},\mathrm{Re}>445\end{array}\right.$$

(31)

$$f=\frac{9}{\sqrt{\mathrm{Re}}}+0.9\frac{0.75E{o}^2}{0.75E{o}^2+4.5}$$

(32)

$$F=\frac{48}{\mathrm{Re}}\left(\frac{1+12M{o}^{1/3}}{1+36M{o}^{1/3}}\right)+0.9\frac{E{o}^{3/2}}{1.4\left(1+30M{o}^{1/6}\right)+E{o}^{3/2}}$$

(33)

$${\left(\frac{a}{R}\right)}^2=\frac{2}{3.974\times {10}^{-3}{\left(We-12.62\right)}^2-7.186\times {10}^{-4}{\left(Eo-17.87\right)}^2}$$

(34)

$${\left(\frac{A}{Ro}\right)}^2=\frac{10\left(1+1.3M{o}^{1/6}\right)+3.1Eo}{10\left(1+1.3M{o}^{1/6}\right)+Eo}$$

(35)

where, $Eo$ is the Eotvos number, dimensionless, $Eo=\frac{g({\rho }_l-{\rho }_g){d_e}^2}{\sigma }$. $d_e$ is the bubble equivalent diameter, m. $\sigma$ is the surface tension of drilling fluid, N/m. $Mo$ is Morton number, dimensionless, $Mo=\frac{g({\rho }_l-{\rho }_g){\mu }^4}{{\rho }_l^2{\sigma }^3}$. $We$ is a Weber number, dimensionless, $We=\frac{({\rho }_l-{\rho }_g)v_g{d}_e}{\sigma }$.

3.4. Model Verification

The clean bubble migration data obtained by Rehder et al. [28] from the observation of methane bubbles released at a depth of 475.9 m were used to verify the clean bubble migration velocity model established in this study. The observation data were collected from 0 to 150 s after the release of bubbles. At the same time, the hydrated bubble migration data obtained by Rehder et al. [28] at the depth of 1051.4 m were used to verify the hydrated bubble migration velocity model. The observed data were collected from 200 to 500 s after the release of bubbles.

Figure 8(a) shows the comparison between the calculated results of the clean bubble migration velocity model and the experimental data. The maximum prediction error of this model is less than 5%. Figure 8(b) shows the comparison between the calculated results of the hydrated bubble migration velocity model and the experimental data. The maximum prediction error of this model is less than 10%. Due to the influence of complex processes such as the dynamic growth of hydrate shell and dissolution mass transfer in the migration process of hydrated bubble, the accuracy of the migration velocity prediction model of hydrated bubble is lower than that of the migration velocity prediction model of clean bubble. The error analysis shows that the bubble migration velocity model established in this study can basically realize the accurate prediction of bubble migration velocity.

4. Discussion

In deepwater wellbore annulus, key parameters such as drag coefficient and bubble equivalent radius in the bubble migration velocity model are affected by the initial bubble size and annulus fluid properties, which also lead to differences in the safe shut-in cycle. Therefore, it is necessary to analyze the influence of initial bubble size, annular fluid viscosity and annular fluid density on the safe shut-in cycle based on the new model of bubble migration velocity. In this study, basic parameters of a deepwater gas well in the eastern South China Sea are used to analyze the influence of different parameters on bubble migration velocity and safe shut-in cycle. The well encountered severe sea conditions when drilling to the reservoir section, and shut-in measures were taken. The main component of reservoir gas is methane. Basic parameters and relevant model parameters are shown in Table 3.

4.1. Influence of initial bubble size

Due to the heterogeneity of reservoir porosity, the size of the initially formed bubbles is different, so it is necessary to explore the migration laws with different sizes of bubbles in deepwater wellbore annulus. Based on the new model of bubble migration velocity, the variation law of bubble migration velocity with well depth and safe shut-in cycle of different initial diameters are calculated under the condition of bottom hole fluid viscosity of 10 mPa·s and density of 1250 kg/m³.

Figure 9 shows the variation of the migration velocity of hydrated bubbles and clean bubbles with well depth when the initial bubble diameter is 2 mm, 4 mm and 6 mm. In the process of clean bubble migration, the velocity from bottom hole to wellhead gradually increases, but the increasing rate gradually slows down, which is the result of the comprehensive effect of bubble volume change, dissolution and fluid viscosity change. In order to recycle drilling fluid and simplify the treatment process of returned drilling fluid during deepwater drilling, the gas solubility of drilling fluid is very small. Therefore, in the process of clean bubble migration, the decrease of pressure leads to the increase of bubble volume, while the effect of dissolution has little effect on the change of bubble volume. In addition, the annular temperature gradually decreases and the drilling fluid viscosity increases accordingly from the bottom of the well to the subsea wellhead, which inhibits the increase of the bubble migration velocity to a certain extent. The migration velocity of hydrated bubble is divided into gradually decreasing stage and slowly increasing stage. When bubble initially enters the hydrate formation area, hydrated bubble is gradually formed due to the growth and coating of the hydrate shell. Gas consumption and hydrate shell thickening play a dominant role, which leads to the gradual decline of bubble migration velocity. With the increase of hydrate shell thickness and the renewal of pores inside the hydrate shell, both the inward mass transfer of water and the outward mass transfer of gas are inhibited. The thickness of hydrate shell gradually increases slowly and the dissolution is limited. The decrease of pressure leads to the increase of bubble volume, which gradually plays a dominant role, and the bubble migration velocity also increases accordingly. In the process of migration near the subsea wellhead, the fluid viscosity gradually increases due to low temperature, which also causes the bubble migration velocity increases slowly. Comparing the velocity change of bubbles with different sizes in the migration process, it can be found that the effect of hydrate transition on the migration velocity of large bubbles is more obvious. When the initial diameter of bubbles is 2 mm and 6 mm, the migration velocity of hydrated bubbles is about 0.010 m/s and 0.031 m/s lower than that of clean bubbles at the wellhead, respectively.

Figure 10 shows the migration cycles of hydrated bubbles and clean bubbles when the initial diameter of bubbles is 2 mm, 4 mm and 6 mm. Since the bubble migration to the subsea wellhead can bring severe challenges to subsequent well opening operations, the bubble migration cycle is generally regarded as the safe shut-in cycle [14]. Comparing the migration cycles of bubbles with different sizes, it is found that the safe shut-in cycle decreases significantly with the increase of the initial bubble diameter. When bubble diameter is 2 mm and 6 mm, the safe shut-in cycle is 42.7 h and 13.6 h respectively. Due to the coating and thickening growth of the hydrate shell, the migration velocity of the bubble is significantly reduced after the bubble enter the hydrate formation area. The formation of hydrated bubble can prolong the safe shut-in cycle. When the initial diameter of the bubble is 2 mm and 6 mm, compared with the clean bubble, the safe shut-in cycle of the hydrated bubble can be extended by 5.1 h and 1.6 h respectively. The above analysis shows that both bubble size and hydrate phase transition have significant effects on bubble migration velocity. The accurate calculation of safe shut-in cycle needs to consider the hydrate phase transition characteristics in deepwater wellbore environment and the change of equivalent diameter during bubble migration.

4.2. Influence of annular fluid viscosity

The above research found that the change of annular fluid viscosity during bubble migration will affect its migration velocity. In order to further explore the influence of different fluid viscosity on bubble migration velocity, the variation law of bubble migration velocity with well depth and safe shut-in cycle are calculated under the conditions of initial bubble diameter of 4 mm and bottom hole fluid density of 1250 kg/m³.

Figure 11 shows the variation of the migration velocity of hydrated bubbles and clean bubbles with well depth when the bottom hole fluid viscosity is 10 mPa·s, 20 mPa·s and 30 mPa·s. With the increase of bottom hole fluid viscosity, the migration velocity of clean bubbles decreases significantly and the migration velocity of clean bubbles in high viscosity fluid increases slowly. When the viscosity of bottom hole fluid is 10 mPa·s and 30 mPa·s, the velocity of clean bubble migration to wellhead increases by 0.021 m/s and 0.004 m/s respectively. This is caused by the significant increase of drag coefficient with the increase of fluid viscosity. For hydrated bubbles, the increase of fluid viscosity shortens the migration distance during the thickening growth of hydrate shell. The migration velocity of hydrated bubbles is more sensitive to the change of fluid viscosity. Due to the lower temperature near the subsea wellhead, the fluid viscosity increases significantly and the migration velocity of hydrated bubbles increases more slowly than that of clean bubbles in this section.

Figure 12 shows the variation of bubble migration cycle in annular fluid with different viscosity. Comparing the migration cycles of clean bubbles and hydrated bubbles under different fluid viscosity, the formation of hydrated bubbles is more conducive to prolonging the safe shut-in cycle with the increase of fluid viscosity. When the viscosity of bottom hole fluid is 30 mPa·s, the safe shut-in cycle of hydrated bubble is 16.8 h longer than that of clean bubble. The above analysis shows that the increase of bottom hole fluid viscosity can significantly decrease the bubble migration velocity and prolong the bubble migration cycle, which is more obvious in the thickening growth section of the hydrate shell. Therefore, the safe shut-in cycle can be greatly improved by means of "double liquid plug injection", which is conducive to dealing with complex downhole accidents and harsh sea conditions. Double liquid plug injection means that a segment of high viscosity liquid is injected into bottom hole reservoir section and hydrate shell growth and thickening section respectively. Injecting a segment of fluid with high viscosity at the bottom of the well can effectively reduce the initial bubble migration velocity. At the same time, injecting a segment of high viscosity liquid into the thickening growth section of hydrate shell can effectively prolong the migration time of bubbles in this section and greatly decrease the migration velocity of hydrated bubbles.

4.3. Influence of annulus fluid density

The annulus fluid density will affect the annulus pressure distribution, which will affect the bubble migration velocity. In order to further explore the influence law of different annulus fluid density on bubble migration velocity, the variation law of bubble migration velocity with well depth and safe shut-in cycle are calculated under the conditions of initial bubble diameter of 4 mm and bottom hole fluid viscosity of 10 mPa·s.

Figure 13 shows the variation of the migration velocity of hydrated bubbles and clean bubbles with well depth when the bottom hole fluid density is 1200 kg/m³, 1250 kg/m³ and 1300 kg/m³. The influence of annulus fluid density on the migration velocity of hydrated bubble and clean bubble is consistent. With the increase of the bottom hole fluid density, the migration velocity of hydrated bubble and clean bubble increases slightly. When the fluid density is 1200 kg/m³, the migration velocity of hydrated bubble and clean bubble to the wellhead is 0.036 m/s and 0.052 m/s respectively. While when the fluid density is 1300 kg/m³, the migration velocity of hydrated bubble and clean bubble to the wellhead is 0.038 m/s and 0.055 m/s respectively. This indicates that the effect of fluid density on bubble migration velocity is not significant. This is consistent with the variation of bubble migration cycle in fluid with different density in Figure 14. When the bottom hole fluid density is 1200 kg/m³ and 1300 kg/m³, the migration cycle of hydrated bubbles is about 2.0 h and 1.9 h longer than that of clean bubbles respectively. There is little difference in safe shut-in cycle under different annulus fluid density. The above analysis shows that the migration velocity of hydrated bubbles and clean bubbles is not sensitive to the change of annular fluid density.

5. Conclusions

In this study, a new model of bubble migration velocity in deepwater wellbore annulus is established considering the effects of hydrate phase transition and gas-water bidirectional mass transfer. Combined with the influence of bubble size and fluid properties in deepwater wellbore annulus on key parameters of bubble migration velocity model, the variation rules of bubble migration velocity and safe shut-in cycle under different initial bubble size, fluid viscosity and density are analyzed. The main conclusions are as follows:

(1) The migration velocity of hydrated bubble is divided into gradually decreasing stage and slowly increasing stage. The gas consumption and the thickening of hydrate shell in the gradually decreasing stage play a dominant role, and the increase of bubble volume caused by the decrease of pressure in the slowly increasing stage is the most important factor.

(2) The formation of hydrated bubble can significantly reduce the migration velocity of bubble and effectively prolong the safe shut-in period. The migration cycle of hydrated bubble can be significantly increased by decreasing bubble size and increasing annular fluid viscosity.

(3) The initial size of bubble and the viscosity of annulus fluid are the main factors affecting the migration velocity of bubble, while the density of annulus fluid has little effect on the migration velocity of hydrated bubble and clean bubble.