The physical model of wellbore fluid flow after gas invasion during offshore MPD is shown in Figure 1. The physical process of fluid flow can be described as follows:

(1) Gas flows from the formation to the wellbore during gas invasion;

This paper establishes a pressure wave propagation model for the process of pressure wave propagation formed by the added backpressure in step (4). The main assumptions include:

(1) The fluid in wellbore and wellbore wall are in a thermal equilibrium state;

The control volumes of gas and liquid phases are taken out separately and then mass conservation equations are established, as shown in Eq. 1. ∂ A α g ρ g ∂ t + ∂ A α g ρ g u g ∂ z = q g − Γ g ∂ A α l ρ l ∂ t + ∂ A α l ρ l u l ∂ z = Γ g (1) where A means area of annulus, m²; α means volume fraction of fluid, dimensionless; ρ is the density of fluid, kg/m³; u means the velocity of fluid, m/s; q g means gas influx rate, k g / m 3 ∙ s ; Γ g means the gas-liquid inter-phase mass transfer rate, k g / m 3 ∙ s ; the subscripts g and l represent gas and liquid, respectively.

The momentum conservation equations for gas and liquid phases are as follows: ∂ A α g ρ g u g ∂ t + ∂ A α g ρ g u g 2 ∂ z + ∂ A α g P ∂ z = F g l ∂ A α l ρ l u l ∂ t + ∂ A α l ρ l u l 2 ∂ z + ∂ A α l P ∂ z = F lg (2) where P is the annulus pressure, Pa; F represents the external force acting on the fluid.

(1) Gas influx rate

During the MPD process, once the bottom hole pressure is lower than the reservoir pressure, the reservoir gas would flow into the annulus in a negative pressure manner. The flow of reservoir gas into the annulus can be treated as a plane radial flow of an infinite formation flowing towards a point, which can be described by the following model, as shown in Eq. 3 (NygaardVefringFjelde et al., 2007): P k = P e − q μ g 4 π k z 2 s + ⁡ ln 4 k t e γ ϕ μ c r w 2 (3) where P _k is the bottomhole pressure when gas influx occurs, Pa; P _e means the reservoir pressure, Pa; μ g is the viscosity of gas, Pa ∙ s ; k is the permeability of reservoir, mD; z means the length of open hole section, m; s means the skin factor of reservoir, dimensionless; t indicates the time of gas influx, s; γ is the Euler’s constant, equals to 0.557721, dimensionless; ϕ is the porosity of the reservoir, dimensionless; r _w is the radius of wellbore, m.

The open hole section of wellbore can be divided into n part, the length of every part is ∆ z , then gas influx rate in each control volume i can be described as: q g i = 4 π k i z i P e − P k μ 2 s + ⁡ ln 4 k i t e γ ϕ μ c r w 2 (4)

(2) Gas-liquid inter-phase mass transfer rate induced by gas dissolution Γ g

Gas density is calculated with the real gas state equation, as shown in Eq. 7. ρ g = M g P R T Z (7) where M g means molar mass of gas, kg/mol, in this paper, it is assumed that all invading gas is methane, therefore M g = 0.016 k g / m o l ; R represents the gas constant, Z represents the gas deviation factor, and Z is calculated using the following equation (Eq. 8): Z = 1 + A 1 + A 2 T pr + A 3 T pr 3 ρ pr + A 4 + A 5 T pr ρ pr 2 + A 5 A 6 ρ pr 5 / T pr + A 7 ρ pr 2 T pr 3 exp − A 8 ρ pr 2 (8) where A₁∼A₈ are the coefficients; ρ pr indicates dimensionless contrast density, T pr represents apparent contrast temperature.

The gas viscosity can be calculated using the following equations (Eq. 9–Eq. 12). μ g = 10 − 4 K e x p X ρ g Y (9) K = 2.6832 × 10 − 2 470 + M T 1.5 116.1111 + 10.5556 M + T (10) X = 0.01009 350 + 54777.7 T + M (11) Y = 2.447 − 0.2224 X (12) where μ g indicates the viscosity of methane at temperature T and pressure P.

(4) Gas-liquid interfacial resistance

The external force acting on the gas phase is shown in Eq. 13. F g l = − M l i d + M l i n d + α g ρ g g ⁡ sin ⁡ θ + 4 τ g w D (13)

The external force acting on the liquid phase is shown in Eq. 14. F lg = M l i d + M l i n d + α l ρ l g ⁡ sin ⁡ θ + 4 τ l w D (14) where M l i d + M l i n d means the gas-liquid interfacial resistance; α g ρ g g ⁡ sin ⁡ θ indicates the gas phase gravity; 4 τ g w D is the shear stress between gas and casing wall or wellbore wall; D represents the equivalent hydraulic diameter of the annulus; M l i d represents the momentum exchange term caused by drag stress, which is expressed as Eq. 15. M l i d = 3 8 C D α g ρ l u r u r r b (15) where C D represents the drag force coefficient, dimensionless; r b represents the average radius of the bubble, m; u r represents the difference in gas-liquid velocity, m/s, expressed as Eq. 16. u r = u g − u l (16) M l i n d represents the momentum exchange term caused by non dragging stress, and its expression is shown in Eq. 17. M l i n d = C v m α g ρ l ∂ u g ∂ t + u g ∂ u g ∂ z − ∂ u l ∂ t + u l ∂ u l ∂ z (17) where, C v m is the virtual mass force coefficient, which is closely related to the void fraction and flow pattern of gas-liquid two-phase flow, the empirical expression for bubbly flow patterns is as follows: C v m = 1 + 2 α g 2 − 2 α g (18)

The mass conservation equations, momentum conservation equations, and auxiliary equations form a gas-liquid two-phase dual fluid model. Based on this model, the propagation velocity of pressure wave in gas-liquid two-phase flow is derived. ∂ α g ρ g ∂ t + ∂ α g ρ g u g ∂ z = 4 π k i z i P e − P k μ 2 s + ⁡ ln 4 k i t e γ ϕ μ c r w 2 − 0.1781 × α l S G g P × 10 0.0125 A P I 18 × 10 0.00091 1.8 T − 459.67 1.20482 ∂ α l ρ l ∂ t + ∂ α l ρ l u l ∂ z = 0.1781 × α l S G g P × 10 0.0125 A P I 18 × 10 0.00091 1.8 T − 459.67 1.20482 ∂ A α g ρ g u g ∂ t + ∂ A α g ρ g u g 2 ∂ z + ∂ A α g P ∂ z = − 3 8 C D α g ρ l u r u r r b + C v m α g ρ l ∂ u g ∂ t + u g ∂ u g ∂ z − ∂ u l ∂ t + u l ∂ u l ∂ z + α g ρ g g ⁡ sin ⁡ θ + 4 τ g w D ∂ A α l ρ l u l ∂ t + ∂ A α l ρ l u l 2 ∂ z + ∂ A α l P ∂ z = 3 8 C D α g ρ l u r u r r b + C v m α g ρ l ∂ u g ∂ t + u g ∂ u g ∂ z − ∂ u l ∂ t + u l ∂ u l ∂ z + α l ρ l g ⁡ sin ⁡ θ + 4 τ l w D (19)

The partial differential components in the above equations are expanded into partial differential forms with V= ( α g , u g , u l , P) as variables, and these equations are organized into the following form. U t ∂ V ∂ t + U z ∂ V ∂ z = W t (20) U t = ρ g 0 0 α g ∂ ρ g ∂ p − ρ l 0 0 α l ∂ ρ l ∂ p ρ g u g α g ρ g + C v m ρ l − C v m α g ρ l α g u g ∂ ρ g ∂ p − ρ l u l − C v m α g ρ l ρ l α l + C v m α g α l u l ∂ ρ l ∂ p (21) U z = ρ g u g α g ρ g 0 α g u g ∂ ρ g ∂ p − ρ l u l 0 α l ρ l α l u l ∂ ρ l ∂ p ρ g u g 2 α g u g 2 ρ g + C v m ρ l − C v m α g ρ l u l α g u g 2 ∂ ρ g ∂ p + 1 − ρ l u l 2 − C v m α g ρ l u g ρ l u l 2 α l + C v m α g α l u l 2 ∂ ρ l ∂ p + 1 (22) W T = 4 π k i z i P e − P k μ 2 s + ⁡ ln 4 k i t e γ ϕ μ c r w 2 − 0.1781 × α l S G g P × 10 0.0125 A P I 18 × 10 0.00091 1.8 T − 459.67 1.20482 − 0.1781 × α l S G g P × 10 0.0125 A P I 18 × 10 0.00091 1.8 T − 459.67 1.20482 3 8 C D α g ρ l u r u r r b + α g ρ g g ⁡ sin ⁡ θ + 4 τ g w D 3 8 C D α g ρ l u r u r r b + α l ρ l g ⁡ sin ⁡ θ + 4 τ l w D (23)

The eigenvalues of the coefficient matrix ( λ ) of Eq. 20 are related to the propagation velocity of pressure waves (Li et al., 2016). Thus, by solving the eigenvalues of the coefficient matrix of the equation system, the propagation velocity of pressure waves in gas-liquid two-phase flow considering gas-liquid mass transfer can be obtained.

In gas-liquid two-phase flow, both gas and water are considered as continuous media. The propagation velocity of pressure waves in each phase of the fluid can be described by the following equation: c = d p d ρ (24)

The pressure wave velocity of two-phase flow considering gas-liquid mass transfer can be obtained by solving the following determinant: ρ g u g − c α g ρ g 0 α g 1 c g u g − c ρ l − u l + c 0 α l ρ l α l 1 c l u l − c ρ g u g u g − c α g u g 2 ρ g + C v m ρ l − c ρ g + C v m ρ l C v m α g ρ l − u l + c α g u g 2 1 c g + 1 − c u g 1 c g ρ l u l − u l + c C v m α g ρ l − u g + c ρ l u l 2 α l + C v m α g − c α l + C v m α g α l u l 2 1 c l + 1 − c u l 1 c l = 0 (25)

The steps for solving the propagation velocity of pressure wave at various positions in the wellbore in gas-liquid two-phase flow are as follows:

(1) Input basic parameters, boundary conditions, and initial conditions;

The reliability of the pressure wave propagation velocity model established in this paper is verified through experimental data. It is very difficult to directly measure the pressure wave velocity transmitted in the gas-liquid two-phase flow of the wellbore during drilling operations. Therefore, the pressure wave velocity measured in vertical gas-liquid two-phase pipe flow experiments by Henry et al. was used to validate the model. Henry et al. used a 2-inch diameter vertical stainless steel tube in their experiment to test the pressure wave velocity in a stable bubble flow with a gas volume fraction less than 40% (Henry et al., 1971).

As shown in Figure 2, the comparison results show that the error between the calculated pressure wave propagation velocity of model in this paper and the experimental data is within 15% with 6% average error. Thus, the accuracy of the pressure wave velocity model in this article is considered acceptable for wave velocity calculation in drilling engineering.

In this subsection, the propagation characteristics of pressure waves at different stages of gas invasion are investigated. Oil based drilling fluid with a density of 1.13 g/cm³ was used, and the difference between bottom hole pressure and formation pressure is set to 0.3 MPa. The invading gas in this paper is 100% pure methane with a density of 0.64 g/cm³. To thoroughly investigate the propagation characteristics of pressure waves in different states of gas invasion, the pit gains of 1 m³, 2 m³, 5 m³ and 10 m³ in this section. It must be noted that this assumed degree of gas invasion in this study is significant. In fact, if the a degree of gas invasion is that significant during the drilling process, this well would be shut down instead of continuing to use the MPD technology. The pressure wave velocity distributions at pit gain of 1 m³, 2 m³, 5 m³ and 10 m³ are calculated and compared.

Figure 3 shows the comparison of pressure wave velocity distributions when pit gains are 1 m³, 2 m³, 5 m³ and 10 m³. As described in Figure 3, the distributions of pressure wave propagation velocity along the wellbore are similar when the pit gain are is 1 m³, 2 m³, and 5 m³. The propagation velocity under these three gas invasion stages above 4,749 m are the same and relatively high. Below 4,749 m, propagation velocities are decreased by the increment of the pit gain. The distribution of pressure wave propagation velocity with a pit gain of 10 m³ is significantly different from these three stages. The low pressure wave propagation velocity length of 10 m³ pit gain is 1,135 m which is highly longer than 50 m of 1 m³, 2 m³, and 5 m³. The times required for pressure wave to propagate from the wellhead to the bottom of the well at these four gas invasion stages has been calculated. Results of times required are 4.3 s, 4.3 s, 4.4 s and 21.8 s, respectively. It can be seen that the time required for pressure waves to propagate in the wellbore changes relatively little under 1 m³, 2 m³, and 5 m³ pit gain. However, the propagation of pressure waves in the wellbore takes a longer propagation time under 10 m³ pit gain. The pressure wave propagates to the bottom of the well to restrain gas invasion, which cannot be completed instantaneously. This indicates that the propagation delay effect of pressure wave must be considered in wellhead back pressure design during well killing stage.

Figure 4 shows the comparison of free gas volume fraction at different gas invasion stage. As depicted in Figure 4, the free gas volume fraction near the bottom of the well increases significantly in the later stage of gas invasion, which is the main reason for the decrease in pressure wave propagation velocity. Due to the ability of oil-based drilling fluids to dissolve gas invading the wellbore, the propagation velocity of pressure waves is less affected by the invading gas under 1 m³, 2 m³, and 5 m³ pit gain. When the gas invasion is small, the gas can dissolve into oil-based drilling fluid through gas-liquid mass transfer. The free gas content in the wellbore is small, and the propagation process of pressure waves in the wellbore is mostly through single-phase liquid. The propagation velocity of pressure waves is relatively high, and the required time is short. As the amount of gas invasion increases, such as the mud pit increment reaches 10 m³, the gas-liquid mass transfer process of the drilling fluid reaches equilibrium, and the dissolution of the drilling fluid is saturated. Therefore, the gas in the drilling fluid becomes more and more abundant. The mass transfer between the gas and liquid phases in the wellbore becomes increasingly strong until it reaches stability. Therefore, an increasing number of gas exists in the form of free gas and reduces the propagation speed of pressure waves.

According to the calculation results and analysis in this section, the control strategy of wellbore pressure in MPD technology needs to be continuously adjusted according to the gas invasion state. The smaller the gas invasion, the earlier it is discovered, and the faster the response of wellbore pressure. The response time of wellbore pressure in the later stage of invasion needs to be calculated based on the gas invasion situation.

In this subsection, the influence of drilling fluid type on pressure wave velocity is investigated. The propagation velocities of pressure wave at different positions in the wellbore under two different drilling fluid are calculated and compared. Figure 5 shows the comparison of pressure wave velocities between water based drilling fluid and oil based drilling fluid when the pit gain is 1 m³. The oil-water ratio of oil-based drilling fluid used in this part is 100:0, while the oil-water ratio of water-based drilling is 0:100.

As shown in Figure 5, the propagation velocity of pressure wave is unevenly distributed at various positions in the wellbore. The propagation velocity of pressure wave in the upper part of the wellbore is relatively high, while the propagation velocity of pressure wave in the lower part of the wellbore is relatively low. In the case of oil-based drilling fluid, from the bottom of well to wellhead, the propagation velocity of pressure wave first decreases from 124.8 m/s to 95.1 m/s, then increases to 113.6 m/s, and then remains at 1,190 m/s. While in the case of water-based drilling fluid, from the bottom of well to wellhead, the propagation velocities of pressure waves increase from 41.7 m/s to 69.5 m/s, and then decrease to 44.4 m/s, and then remains at 1,190 m/s. As indicated by Figure 5, the location of the minimum propagation velocity of wellbore pressure wave in the case of oil-based drilling fluid and water-based drilling fluid is 20 m and 1,135 m away from the bottom of the well respectively.

The influence of drilling fluid type on the propagation velocity of pressure wave is due to the influence of the existing form of invading gas in the wellbore. As shown in Figure 6, in gas-liquid two-phase fluid, the increase of free gas volume fraction would reduce the propagation velocity of pressure wave. In oil-based drilling fluids, most of the invading gases exist in the drilling fluid in liquid form due to the mass transfer between gas and liquid phases. On the contrary, in water-based drilling fluid, there is no inter-phase mass transfer behavior between the invading gas and the drilling fluid. The invading gas exists in the form of free gas, which greatly reduces the propagation speed of pressure waves in gas-liquid two-phase fluid. The distribution of free gas volume fraction in the annulus under two types of drilling fluids is calculated and shown in Figure 6.

Figure 6 shows the comparison of free gas volume fractions in annulus between oil-based and water-based drilling fluids when pit gain is 1 m³. As shown in Figure 6, the free gas volume fraction in the annulus of oil-based drilling fluid is basically close to zero, while the free gas volume fraction in water-based drilling fluid is high. The free gas volume fraction in the annulus of water-based drilling fluid is unevenly distributed. From the bottom to the wellhead, the free gas volume fraction increases from 1% to 11.4%, and then decreases to 0%. The maximum free gas volume fraction is located 78.5 m away from the bottom of the well. The form of invading gas in the annulus has an impact on the density of the annulus mixture. To analyze the effect of annular mixture density on pressure wave propagation velocity, the dissolved gas mass fraction and mud density distribution of annular mixture under two types of drilling fluids are calculated and displayed in Figure 7.

Figure 7 shows the comparison of dissolved gas mass fraction mud density in annulus between oil-based and water-based drilling fluids when the pit gain is 1 m³. As depicted in Figure 7, the mass fraction of dissolved gas in the annulus is 0 in the case of water-based drilling fluid, while in the case of oil-based drilling fluid, the dissolved gas mass fraction in the annulus are relatively high. From the bottom of the well, the dissolved gas mass fraction increases from 0.16% to 2.1%, and then decreases to 0. The maximum position of the dissolved gas mass fraction is 50.4 m from the bottom of the well. From the bottom of the well, the density of oil-based drilling fluid decrease from 1.130 g/cm³ to 1.075 g/cm³, then increased to 1.2 g/cm³, then decreased upwards to 1.115 g/cm³, and finally increased to 1.118 g/cm³. Combining Figures 7A, B, dissolved gas mass fraction reduces the density of oil-based drilling fluid. Based on the comprehensive analysis of Figures 5–7, it can be concluded that the influence of drilling fluid type on the propagation velocity of pressure wave is mainly controlled by the free gas volume fraction, while the influence of annular mud density is relatively small. The response time to wellbore pressure after water based drilling fluid gas invasion is longer. This makes it more difficult to regulate the pressure of water-based drilling fluid after gas invasion than oil-based drilling fluid.

In this subsection, the influence of drilling fluid density on the propagation velocity of pressure waves is studied. The drilling fluid type is oil-based with different fluid density. The propagation velocities of pressure waves in annulus at 200 s of gas invasion under three different drilling fluid densities are calculated and compared.

Figure 8 shows the comparison of pressure wave propagation velocity distribution after 200 s of gas invasion at drilling fluid densities of 1.11 g/cm³, 1.13 g/cm³, and 1.15 g/cm³. As depicted in Figure 8, after gas invasion for 200 s the pressure wave velocity above the well depth of 4,731 m under three different drilling fluid densities are greater than that below 4,731 m. Moreover, the propagation velocities of pressure wave below 4,731 m at a density of 1.15 g/cm³ is greater than that at densities of 1.11 g/cm³ and 1.13 g/cm³. To illustrate the mechanism of the influence of drilling fluid density on the propagation velocity of pressure wave, the distribution of free gas volume fraction in the annulus under three different drilling fluid densities is calculated and displayed in Figure 9.

As shown in Figure 9, the distribution characteristics of free gas volume fraction in the annulus vary under different drilling fluid densities. The free gas volume fraction is unevenly distributed along the wellbore under three different drilling fluid densities. There is a peak free gas volume fraction at the bottom of the wellbore. The increase in drilling fluid density can reduce the peak value of free gas volume fraction at the bottom of the wellbore. This is mainly due to the fact that high drilling fluid density can suppress gas invasion and reduce the volume of gas invading the wellbore. In addition, high drilling fluid density can also enhance mass transfer between gas-liquid phases, thereby reducing the volume fraction of free gas. It can be concluded that under high drilling fluid density, the response and control time for wellhead backpressure are shorter than those of low drilling fluid density.

In this subsection, the propagation velocity of wellbore pressure wave under different backpressure values are calculated and compared to investigate the Influence of backpressure value on the pressure wave velocity.

Figure 10 shows the distribution of pressure wave propagation velocity at 300 s, 600 s, 900 s, and 1,200 s after the occurrence of gas invasion under 0.2 MPa and 0.5 MPa backpressure. As depicted in Figure 10, when a gas invasion occurs for 300 s, there is little difference in the distribution of pressure wave propagation velocity between 0.2 MPa and 0.5 MPa backpressure. In both cases, the propagation velocity of pressure wave near the bottom of the well is low due to the presence of free gas, while the propagation velocity of pressure wave at the upper position is high. Comparing the four figures in Figure 10, the length of the low pressure wave propagation velocity region near the wellbore increases over time. This is because gas invaded into wellbore increases with time, and it is transported from wellbore to wellhead with circulated drilling fluid. Compared with 0.2 MPa backpressure case, the wellbore length with lower pressure wave propagation velocity under 0.5 MPa backpressure is significantly shorter. The reason behind this phenomena is that high backpressure can effectively restrain the velocity of gas invasion into the wellbore, and therefore reduce free gas volume fraction in the annulus.

Figure 11 shows the comparison of time required for pressure waves to propagate from the wellhead to the bottom of the well under different backpressure values. As depicted in Figure 11, the propagation time of pressure waves under 0.2 MPa is greater than that under 0.5 MPa, and the difference between these two cases increase with the increase of gas invasion time. The difference between these two cases increases from 1.38 s at gas invasion 300 s–37.93 s at gas invasion 1,200 s. This also indicates that the propagation time of pressure waves should be considered in the design of wellhead backpressure values in MPD technology. The propagation time of pressure wave at 1,200 s under 0.5 MPa backpressure is lower than that at 900 s. This is because that the gas invasion under 0.5 MPa is controlled properly and the gas invasin is stopped. At the same time, a part of the invasion gas has been transported out at 1,200 s. Therefore, the amount of free gas at 1,200 s is smaller than that at 900 s.