Numerical investigation of gas‐liquid displacement between borehole and gassy fracture using response surface methodology

Gas‐liquid displacement occurs often in fractured gas reservoirs, and can cause gas kick and mud leakage, resulting in a very high risk of losing well control. To analyze gas‐liquid displacement between borehole and gassy fracture, we used computational fluid dynamics to simulate its behaviors. We also used response surface methodology (RSM) to design numerical experiments. The effects of fracture width, bottom‐hole differential pressure, mud density, mud viscosity, and mud displacement were taken into account. We used RSM to determine the influence of the multifactor interaction of gas‐liquid displacement and established an empirical formula for the gas displacement rate. The results show that gas‐liquid displacement is proportional to fracture width, bottom‐hole differential pressure, mud density, and mud displacement; however, the displacement is inversely proportional to mud viscosity. The sensitivity sequence of the gas‐liquid displacement rate is fracture width > bottom‐hole differential pressure > mud viscosity > mud density > mud velocity. The impact of fracture width is clearly higher than that of the other factors, while the mud velocity has almost no impact. Our established empirical formula can be used to predict bottom‐hole gas kick and drilling mud leakage and to inversely predict the fracture width and formation gas pressure.


| INTRODUCTION
As shallow oil and gas reservoir resources are gradually depleted, the inevitable trend for drilling operations is to drill into much deeper formation, among which carbonate formation is the most representative formation. In China, marine carbonate rock is rich in oil and gas. The predicted petroleum geological reserves are approximately 340 × 10 8 t, natural gas geological reserves are approximately 24.3 × 10 12 m 3 , accumulated proved petroleum geological reserves are approximately 29.34 × 10 8 t, and natural gas geological reserves are approximately 3.37 × 10 12 m 3 . 1 Fractures and fracture cavities develop in carbonate formations. Because of the different densities of mud (or drilling fluid) and natural gas, significant gas-liquid gravity displacement occurs in the fracture, and it is a major safety problem for drilling operation in carbonate | 741 formations. It can lead to kick, which if not controlled results in a blowout accident and very large economic losses: a significant challenge to drilling operations and well control. 2,3 Therefore, it is important to study the law of gas-liquid gravity displacement, which can provide scientific guidance for rapid and efficient drilling, prevention of kick and lost circulation.
Much research has been done on the law of mud loss and kicks in fractured formations. Dyke et al 4 proposed a mudloss equation using Newtonian fluid. Their equation could quantitatively analyze mud loss and distinguish the matrix permeability or the fractures by fluid losses. Sanfillippo et al 5 presented a model based on a pressure diffusivity equation that could be used to calculate the fracture width of a single conductive fracture. Lietard et al 6,7 developed a model based on the radial flow of a Bingham plastic fluid into an unlimited extension fracture. They believed that the mud losses would eventually stop because of the yield stress of drilling fluid. Majidi et al [8][9][10] presented a mathematical model for Herschel-Bulkley drilling-fluid losses in naturally fractured formations, the effect of mud rheological properties was investigated, and it was found that yield stress can control the ultimate volume of losses, and the shear-thinning effect can greatly decrease the rate of losses.
Ozdemirtas et al 11 studied the effects of fractures in natural geologic formations on the rate of fluid loss and gain rate between the borehole and the fractured formation. The study found that increased fracture surface roughness and partially closing the walls of the fracture reduced the mud gain and loss flow rates. Li et al 12 constructed a transient mathematical model that combined kick with mud loss in a single fracture. The model indicated that the loss rate from pressure at the casing shoe initially increases linearly over time and then decreases exponentially, ultimately becoming constant. Numerical simulations of gravity gas-liquid displacement have been done using computational fluid dynamics (CFD). [13][14][15] The law of gravity displacement under the conditions of shut-in and in normal drilling conditions has also been analyzed. 12 Fracture width, differential pressure, and mud properties are the major factors affecting fluid velocity in a single fracture.
The Rogaland Research Kick Simulator (RF-Kick) has been introduced to simulate lost circulation, and the behavior and reaction to killing operations were supervised during a kick situation. 16 Some laboratory experiments of gas-liquid gravity displacement have been done with similar experimental devices. [13][14][15][17][18][19][20] The devices consisted of three parts: a wellbore-fracture-formation system, a circulation system, and a data acquisition system. The wellbore-fracture-formation system included a wellbore module, a fracture module, and a formation module. Many scholars believe that the fracture width is more sensitive to the gas displacement rate and volume and the mud viscosity clearly inhibits gas-liquid gravity displacement. 16,18,19 The relationships between mud density, fracture width, and mud viscosity with displacement amount are linear, logarithmic, and exponential, respectively. 20 We have found that foreign scholars focus on studying drilling-fluid loss and kick in wellbores without considering the conditions of gas leakage and kick occurring simultaneously in a fracture. Many experiments have qualitatively determined the gas-liquid displacement law by experimental results, but there has been no quantitative analysis of the gas-liquid displacement law. Single-factor analysis had been conducted by experiments, but the impact of multifactor interaction has not been analyzed. Therefore, our study focused on the gas-liquid gravity displacement law and the impact of multifactor interaction.
Computational fluid dynamics numerical simulation for experimental condition is very laborious. In order to study the impact of multifactor interaction, experimental design method was developed for time-saving and workload reduction. There are so many experiment design methods, such as the response surface methodology (RSM), the uniform experiential design, the orthogonal experimental design, and the factorial experimental design. [21][22][23][24] In this study, we used RSM to study the gas-liquid displacement law. The RSM can obtain the effect of multifactor interactions and reduce the number of experiments. Besides, the RSM has highly accurate regression equations with good prediction performance. This paper is organized as follows: In Section 2, we describe the causes and conditions of gravity displacement, establish a set of governing equations for gravity displacement, and propose a computational method that combines boundary conditions and a finite element mode. We also validate the numerical simulation experiment. In Section 3, we describe the principle and characteristics of RSM, we consider the bottom-hole differential pressures, mud density, mud viscosities, mud speeds, and fracture widths, and the numerical experiment was designed and conducted by using RSM. In Section 4, we analyze the CFD simulation results, propose an empirical equation of the gas displacement rate, and analyze and discuss the law of gas displacement volume and rate.

| Causes and conditions of gravity displacement
There are three conditions under which gas-liquid displacement occurs in fractured gas reservoirs 17 : (a) The formation has fracture passages that allow mud to flow from the wellbore into the formation. (b) There is sufficient space in the formation to accommodate the mud loss. (c) The wellbore pressure is located in the gravity displacement window. As shown in Figure 1, the physical model of the gas-liquid displacement consists of drill string, annulus, concentric wellbore, and vertical fracture. According to the pressure balance relationship between the wellbore and formation, three kinds of potential drilling conditions are considered 20 as follows: • Gas kick: • Coexisting gas kick and mud leakage: where p w1 is the wellbore pressure at the top of the fracture, Pa; p w2 is the wellbore pressure at the bottom of the fracture, Pa; and p fg is the gas pressure in the fracture, Pa.

| CFD mathematical models
Based on the properties of the formation fluid and the geometrical characteristic of the fracture, the assumptions were given as follows: (a) Gas is compressible, while the compressibility of liquid is ignored; (b) no mass transfer between liquid and gas; and (c) the wellbore passes through the vertical fracture. A coordinate system is established in the proposed mathematical model, as shown in Figure 1, where the x-axis is the direction of fracture extension, the y-axis is perpendicular to the fracture surface, and the z-axis is along the vertical direction.
The gas-liquid two-phase flow models include the fracture flow and annular flow. Both of these two kinds of flow models consist of a continuity equation and a momentum equation. The gas-liquid interface is determined by the volume-of-fluid (VOF) method. Therefore, the volume concentrations α l and α g are introduced to define the liquid region and gas region, respectively. The volume fraction continuous equation is solved to track the interface between phases.

The mass conservation equation
The gas-fluid two-phase flow must conform to the mass conservation equation, which can be expressed as 25 follows: where α l and α g are given as follows: where ρ i is the density of phase i, kg/m 3 ; u i is the velocity of phase i in the x direction, m/s; v i is the velocity of phase i in the y direction, m/s; w i is the velocity of phase i in the z direction, m/s; t is time, s; Q i is the flow rate of phase i, m 3 /s; α i is the volume concentration of phase i, dimensionless; subscript i = l represents fluid; subscript i = g represents gas; V e is the effective volume of a single unit, m 3 ; and V g is the gas volume inside a unit, m 3 .
In general, Q l = 0 and Q g ≠ 0 mean drilling into gas productive layer, while Q l = 0 and Q g = 0 mean drilling into nonproductive layer.
As the gas is compressible, the gas volume inside a unit can be determined using the gas state equation: where R is gas constant, R = 8.314 J/(mol K); T is temperature, K; Z is gas compressibility factor, dimensionless; V g is the gas volume inside a unit, m 3 ; n is the gas amount of substance, mol; and p g is the gas pressure, Pa.

The momentum conservation equation
Momentum conservation equations are the significant for both fluid and gas, when gas-liquid two-phase flow in the wellbore, the component of gravity in both x and y directions can be ignored, and then, the momentum conservation equations can be expressed as follows 26-28 : where p i is the pressure of phase i, Pa; μ i is the viscosity of phase i, mPa s; g is the gravitational acceleration, m/s 2 ; subscript i = l represents fluid; and subscript i = g represents gas.

| Mathematical model for fracture flow
Due to the small width of fracture, the gas-liquid two-phase flow in the fracture can be simplified as a parallel plate flow. Thus, the velocity in y direction can be ignored for gas-liquid two-phase flow in the fracture, while velocity in both x and z directions cannot be ignored. The component of gravity in both x and y directions can also be ignored, but the influence of boundary-layer effect on gas-liquid two-phase flow in fracture cannot be ignored. Therefore, the continuity equation and momentum equation can be simplified as follows: 1. The mass conservation equation The gas-fluid two-phase flow in fracture still must conform to the mass conservation equation, due to v i = 0, and thus, the mass conservation equation can be expressed as follows:

The momentum conservation equation
Due to v i = 0 and F x = F y = 0, the momentum conservation equation can be expressed as follows: where τ yx , τ zx , τ xz , and τ yz are given as follows: where τ yx , τ zx , τ xz , and τ yz represent the components of shear stress, Pa.
The governing equation has nonlinear characteristics so that the CFD software was used to solve the governing equation. In order to discretize the fluid region to solve the flow field, the PRESTO (pressure staggering option) format was adopted to discretize the pressure terms, and the discretization of the momentum equations was implemented using a second-order upwind difference scheme. Meanwhile, the PISO (pressure implicit with splitting of operators) algorithm was adopted to solve the coupling equations of pressure and velocity equation. Besides, the convergent criteria for all calculations were set as that the residual in the control volume for each equation is less than 10 −5 .

| CFD modeling
We established the physical model of gas-liquid gravity displacement based on the theory of fluid mechanics and finite element analysis. [29][30][31][32] As shown in Figure 2, the diameter of the wellbore was 152 mm, and the diameter of drill string was 88.9 mm. Based on calculation requirements and the actual formation, we set the fracture width, length, and height as 5 mm, 12 m, and 5 m, respectively. The hexahedral mesh grids were used for computational domain, which is able to increase accuracy in resolving the boundary-layer flow. The number of meshes was 1.09 million. The minimum mesh Finite element calculation of model and mesh volume was 4 × 10 −7 m 3 , and the maximum mesh volume was 6 × 10 −6 m 3 . The cell dimensions of the first layer were between 0.5 mm and 1 mm in each direction, and the growth factor was 1.2. The VOF method and the PISO algorithm were used to solve the governing equations. The laminar flow was considered in the fracture, and the gas was monitored in the wellbore. The drilling mud and methane gas were used as fluid in the cases; the mud density ranged 1.07-1.27 g/cm 3 ; the mud viscosity ranged 15-35 mPa s; the mud velocity ranged 1.0-3.0 m/s; and the fracture width ranged 5-15 mm.
Boundary conditions used in the calculation were as follows: 1. Only the mud flowed into the inlet. The inlet-velocity boundary was used for the import of the simulated model, in which the velocity varied from 1.0 to 3.0 m/s. 2. The pressure outlet boundary was used for the export of the simulated model, in which the pressures varied from 29.5 to 30.5 MPa. 3. A no-slip boundary and a standard wall function were used for the model wall. 4. In the initial state, the fractures were filled with methane gas and the pressure was 29 MPa, and the wellbore was filled with mud.

| Model validation
In order to demonstrate the performance of the CFD modeling, the experimental results of gas-liquid gravity displacement published by Shu 17 were utilized. The schematic diagram of experimental device is shown in Figure 3, and the experimental device consists of wellbore, drill string, annulus, fracture, and air supply device, where the wellbore and air supply device were connected by a vertical fracture. The wellbore has a height of 2 m and an inner diameter of 140 mm; the drill string has a height of 2.5 m, an outer diameter of 63 mm, and an inner diameter of 57 mm; the fracture has a length of 600 mm, a height of 400 mm, and a fracture width of 0.5 mm; the air supply device has a height of 0.55 m, an outer diameter of 65 mm, and an inner diameter of 55 mm. Water and air were used for conducting the gas-liquid gravity displacement experiment. The carboxyl methylcellulose (CMC) was used to change water viscosity, the liquid density is 1.0 g/ cm 3 , and the liquid viscosity is 17 mPa s. The outlet of wellbore was connected with atmospheric pressure, and the differential pressure control was achieved by air supply devices in which the pressures were set 10-13 kPa. In the initial state, the fracture was filled with air, and the wellbore was filled with liquid. We compared the numerical simulation results with the experimental results, as shown in Figures 4 and 5.
The results show that the simulation results were close to the actual experimental results; the error was generally less than 15%. This shows that the simulation results can be utilized to calculate gas-liquid gravity displacement in fractured gas reservoirs.

DESIGN USING RSM
Many factors influence the gas-liquid displacement law. We used the RSM to reduce the number of experiments and investigate the impact of multifactor interaction on gas-liquid displacement.

| Principle and characteristics of RSM
Response surface methodology is an experimental design method proposed by Box. 33 It is an optimization method for designing experimental and mathematical modeling. The functional relationship is established by experimenting with representative local points, and the optimal level of each factor is obtained. 34 In engineering optimization design, the RSM obtains not only the changing relationship between the target and the variable but also an optimization scheme.
The RSM also has the advantages of reduction in the number of experiments, high accuracy of regression equations, and good prediction performance. Current methods of constructing response surfaces include mainly polynomials, exponential functions, logarithmic function fitting, and approximation methods such as neural networks. According to the Weierstrass polynomial best approximation theorem, many types of functions can be approximated by polynomials. Therefore, we used the polynomial approximation model to analyze the relationships between variables and functions. This model has the characteristics of simple mathematical expressions, small calculations, and fast convergence, and can be expressed explicitly. It is the most widely used method. Low-order polynomial approximations are usually used within a range of design variables, such as linear or second-order models 35,36 : where Y is the predicted response value; x i and x j are variables; β 0 is an offset term; β i is a linear offset; β ii is a second-order offset coefficient; β ij is an interaction coefficient; k is the number of design variables; and e is the error.

| Experiment design and results
Many factors influence gas-liquid displacement, and the displacement is complex. Therefore, the experiments were designed by RSM to reduce the number of experiments. This paper focuses on five design factors that influence gas-liquid displacement, as shown in Table 1.
The experiment was designed by the RSM, and the results are shown in Table 2. There were 46 group representative experiments to study the gas-liquid gravity displacement, and the experimental results were converted into the gas displacement rate at the standard atmospheric pressure. The differential temperature was ignored.

| Influence of fracture width
Varying fracture widths make the fluid flow area vary during the gas-liquid displacement process, which could affect the formation fluid flow. Therefore, we experimented with gasliquid displacement at three fracture widths (1 mm, 5 mm, and 10 mm), as shown in Figure 6. We also analyzed the effects on the gas-liquid displacement law. We found that the gas-liquid gravity displacement volume, the area of the liquid, and the mud-loss circulation increased as the fracture width increased.

| Influence of bottom-hole differential pressure
During the drilling process, the bottom-hole differential pressure is unstable. This may affect the gas-liquid displacement law, so we experimented with three bottom-hole differential pressures to determine their effect on the gasliquid displacement law. The results are shown in Figure 7. We found that the gas-liquid displacement phenomenon was more pronounced at higher bottom-hole differential pressures. Both the area of the liquid in the fracture and the mud leakage increased as the bottom-hole pressure increased.

| Influence of mud density
Mud density in a wellbore varies, which could affect the gasliquid displacement law. Therefore, we experimented on the gas-liquid displacement at three mud densities to determine the effect on the gas-liquid displacement law, and the results are shown in Figure 8. We found that as the mud density increased, the gas-liquid gravity displacement volume, the area of the liquid, and mud leakage increased. However, the analysis result shows that the gas-liquid displacement phenomenon was not significantly altered with the mud density increase.

| Influence of mud velocity
During the drilling process, the mud velocity in the wellbore varies, which could affect the gas-liquid displacement law. Therefore, we experimented with three mud velocities and their effects on the gas-liquid displacement law. The results are shown in Figure 9. As the mud velocity increased, the area of the mud in the fracture, the gas height in the wellbore, and the mud leakage increased. However, the analysis result shows that the influence of the mud velocity on the gas-liquid displacement phenomenon was insignificant.

| Influence of mud viscosity
Mud viscosity in the wellbore varies, which could affect the gas-liquid displacement law. Therefore, we experimented with three mud viscosities to determine the effect on the gas-liquid displacement law. The results are shown in Figure 10. As the mud viscosity decreased, the volume of the gas-liquid gravity displacement, the area of the liquid, and the mud leakage increased. The analysis result shows that the gravity displacement phenomenon is significant with decreased mud viscosity.

| Analysis of RSM
We used the numerical simulation results shown in Table 2 to propose the following fitted regression equation of gas displacement rate based on the RSM theory: where Q is gas flow, m 3 /s; a is the fracture width, mm; P is the bottom-hole pressure, MPa; ρ is the mud density, g/cm 3 ; v is the mud velocity, m/s; and μ is the mud viscosity, mPa s. We analyzed the regression equation. The adjusted R 2 was a modification of the R 2 that was used to adjust the number of explanatory terms in the model. The predicted R 2 indicates the ability of the model to predict results. The higher value of the predicted R 2 indicates a higher prediction accuracy. The results are shown in Table 3. The R 2 is 0.9958 that the model can accurately explain 99.58% of the experiments in the experimental range. The adjusted R 2 and predicted R 2 are 0.9924 and 0.9839, respectively, which shows that the model is accurate and the error is small. The CV is 0.92%, which shows that the numerical simulation experiment had high reliability and accuracy. The P-value was less than .05, which had a significant impact on the model. The fracture width, bottom-hole differential pressure, and mud viscosity had a significant effect on the gas displacement rate, as shown in Table 3.
We analyzed the relationship between the experimental results and the model prediction results, the residual normal distribution, the residual value of the fitted value, and the residual of each group. The results are shown in Figure 11. The predicted results of the model and the experimental results are concentrated on a straight line. The normal probability of the residual is concentrated on a straight line, and the residual distribution of the predicted value is disorderly and irregular. The error analysis shows that the model has good adaptability, and the prediction error of the gas displacement rate is small in the experimental range.

| Analysis of multifactor influence
The gas-liquid gravity displacement of fractured gas reservoirs is usually controlled by a variety of influencing factors. Therefore, it is necessary to study the gas-liquid displacement law of the influencing factors and analyze the effects of fracture width, mud density, mud viscosity, mud velocity, and various bottom-hole pressures. Consequently, we analyzed the interaction of various factors with the gas-liquid gravity displacement law.

| Variation of gas displacement volume
We combined the RSM with the numerical simulation method and used a 10 mm fracture width as the calculation model. The interaction influence of bottom-hole pressure, mud density, mud viscosity, and mud velocity on the gas displacement volume was investigated. The results are shown in Table 4. The gas displacement volume is the largest when the bottom-hole pressure, mud density, and mud velocity are the maximum value. The gas displacement volume was positively correlated with fracture width, bottom-hole differential pressure, mud density, and mud velocity. In other words, the gas-liquid gravity displacement phenomenon was significant when the bottom-hole differential pressure, mud density, and mud velocity increased. The mud viscosity interacted with other factors, and the gas displacement volume was negatively correlated with it. The gas displacement volume was large with a reduced mud viscosity. The analysis results show that bottom-hole differential pressure and mud viscosity have

| Variation of gas displacement rates
The interactions of fracture width, bottom-hole differential pressure, mud density, mud viscosity, and mud velocity on the gas displacement rate were analyzed by RSM. First, we analyzed the relationship of the gas displacement rate with fracture width, bottom-hole pressure, mud density, mud velocity, and mud viscosity. The results are shown in Figure  12. The gas displacement rate increased with increases in the fracture width, bottom-hole differential pressure, mud velocity, and mud density. The gas displacement rate decreased when the mud viscosity increased. The analysis results show that the fracture width had the most significant effect on the gas displacement rate. The bottom-hole differential pressure and mud viscosity had a great influence on the gas displacement rate. Mud velocity and mud density had little effect on the gas displacement rate. We clarified the gas displacement rate law concerning the influence of single factors on the gas displacement rate, but the interactions of various influencing factors on the gas displacement rate were still unclear. Therefore, we determined the results of the interactions of the factors; the results are shown in Table 5.
The results show that the interaction between the fracture widths and other factors has the most obvious effect on the gas displacement rate, and the gas displacement rate changes mostly with the fracture width change. The interaction between the mud density and the mud velocity has little effect on the gas displacement rate. When the mud viscosity interacts with other effect factors, the interaction with the fracture width has a significant effect on the gas displacement rate. The interaction with the bottom-hole differential pressure has a great influence on the gas displacement rate and the interaction with the mud rate. We analyzed the gas displacement rates against multiple factors and found that the law was similar to the influence of single factors. In the other words, the sensitivity sequence of gas displacement rate was fracture width > bottom-hole differential pressure > mud viscosity > mud density > mud velocity, among which the impact of fracture width was obviously higher than the others, while the mud velocity had almost no impact.

| CONCLUSIONS
We described the causes and conditions of gravity gas-liquid displacement between borehole and gassy fracture, and establish a CFD model for gas-liquid displacement. The numerical simulation was validated using indoor experimental result. The RSM was used to design gas-liquid displacement experiments with various factors, which greatly reduced the duration and number of experiments. The influence of fracture width, bottom-hole differential pressure, mud density, mud velocity, and mud viscosity was considered. An empirical equation of the gas displacement rate was proposed by the RSM, and the law of gas displacement volume and rate were also discussed. The following conclusions can be drawn: 1. We analyzed the CFD results and found that the gas-liquid displacement phenomenon is increasingly pronounced with increases in the fracture width, bottom-hole differential pressure, mud density, and mud velocity. The gas-liquid displacement phenomenon was also more pronounced with a decrease in mud viscosity. 2. We analyzed the gas-liquid displacement law of fractured gas reservoirs using RSM. The results show that the bottom-hole differential pressure and mud viscosity have an appreciable impact on gas displacement volume, but the mud density and mud velocity have an insignificant effect on gas displacement volume. We also analyzed the law of gas displacement rate. The results show that the sensitivity sequence of gas displacement rate is fracture width > bottom-hole differential pressure > mud viscosity > mud density > mud velocity. The impact of fracture width is obviously higher than the others, while mud velocity had almost no impact. Therefore, gas kick and mud leakage can be reduced in fracture formations by properly increasing the mud viscosity and reducing the bottom-hole differential pressure. gas displacement rate in the experimental range. We used the model describing the mud invasion into a single fracture to invert the mud-loss data in order to estimate fracture width and formation gas pressure. However, this formula has some limitations, and more research is required.