A Novel Geothermal Wellbore Model Based on the Drift-Flux Approach

1. Introduction

Geothermal energy, characterized by its low carbon footprint and abundant reserves, has emerged as a highly promising renewable energy source, garnering ever-growing interest. However, as drilling depths escalate, the operational environment within geothermal wells deteriorates significantly, posing challenges for accurate parameter measurements. Consequently, Therefore, employing numerical simulation for modeling geothermal wells is an excellent choice. It allows for an in-depth study of the complex subsurface environment, enabling the prediction of parameters such as bottom hole pressure and temperature. This predictive capability is instrumental in optimizing drilling parameters to ensure safe drilling practices.

In recent years, many researchers have conducted modeling of geothermal wells. Bjornsson et al. have developed HOLA, a steady-state simulator for geothermal wells [1]. Tian et al. have developed a geothermal wellbore model based on the two-phase flow conservation equations. This model is used to simulate multiphase flow inside the wellbore and the heat transfer process between the well and the formation. It is capable of handling complex drilling scenarios and achieves this functionality through input data [2]. Because steady-state geothermal wellbore models cannot accurately reflect the actual flow conditions inside the wellbore, many scholars have begun to use transient models to model geothermal wellbores. Khasani et al. used numerical simulation to establish a wellbore model to investigate the presence of CO2 gas in geothermal water and the effect of scale deposition formed in the wellbore on the flow characteristics of fluids in geothermal wells [3]. Gao et al. developed a coupled THM(Thermo-Hydro-Mechanical) semi-numerical model for fractured reservoirs with three types of flow modes based on porous media theory, including non-local thermal equilibrium and non-Darcy law. The study analyzed the impact of CO2 thermal properties on heat transfer and flow inside the wellbore [4]. The drift-flux model is commonly used to describe multiphase flow, particularly in cases of gas-liquid two-phase flow, making it especially significant for modeling flow in geothermal wells. The drift-flux model was first proposed by Zuber and Findlay [5] in 1965. Hasan et al. used the drift-flux method to establish a robust model for two-phase flow in geothermal wells [6]. Akbar et al. developed a model that utilizes the finite element method to solve heat and fluid flow in high enthalpy deep geothermal wellbores. Additionally, they employed the drift-flux model to simulate the transient heat flow of a compressible, two-phase fluid traveling along the wellbore [7]. Tonkin et al. conducted a review of mathematical models for simulating geothermal wellbores, detailing parameter definitions and differences between models. In 2023, they developed a transient geothermal wellbore simulator based on a mathematical model comprising three conservation equations (mass, momentum, and energy), and incorporating the drift-flux model established by Shi et al [8,9,10]. Lei et al. used the drift-flux model to establish a geothermal wellbore model for controlling two-phase flow containing CO2 in the wellbore. They proposed a robust computational method and applied this model to three production wells in the Yangba geothermal field [11]. Xu et al. developed a transient non-isothermal two-phase flow model for dynamic simulation of multiphase flow in gas kick wellbores, with the ability to predict the temperature and pressure distribution inside the wellbore [12]. Hajidavalloo et al. developed conservation equations for the steady-state annulus and drill string, solved using the finite difference method. They analyzed the impact of temperature changes on predicting bottom hole pressure [13].

From the above literature, it is evident that many researchers have modeled wellbores. However, steady-state wellbore models often have low accuracy and cannot adequately reflect the true conditions inside the wellbore. Some geothermal drilling models include the energy conservation equation, which is essential. Although some models can simulate complex multiphase flows, their descriptions of heat transfer behavior within the wellbore are insufficient, leading to less accurate predictions. Due to inaccuracies in temperature predictions, and since temperature can significantly affect the properties of the fluid in the well, pressure predictions also have inaccuracies.

This paper establishes a comprehensive geothermal wellbore model using the conservation equations of mass, momentum, and energy, combined with Shi et al.'s drift-flux model [10] and Hasan et al.'s heat transfer model [14]. The highly nonlinear equations are discretized using the finite difference method and solved using the Newton-Raphson iteration method. The model is analyzed and validated through numerical simulation software, demonstrating promising application prospects. Furthermore, comparisons with models mentioned in the literature highlight the accuracy and superiority of this model.

2. Model formulation

As shown in Figure 1, the drilling fluid exits through the drill bit nozzles and enters the wellbore annulus, carrying bottom-hole cuttings and gas-liquid mixtures back to the surface while exchanging heat with the formation in the annulus. To establish the geothermal wellbore model, the following assumptions are made:

• The gas-liquid mixture flows one-dimensionally along the wellbore;
• Disregard axial heat conduction within the wellbore;
• Disregard the impact of drill cuttings;
• Consider the impact of the casing and cement sheath on heat transfer;
• Consider the slip relationship between the gas and liquid phases.

2.1. Conservation equations for gas-liquid two-phase flow

In geothermal drilling, gas-liquid two-phase flow must satisfy the conservation laws of mass, momentum, and energy. The equations are as follows:

Mass conservation equation

$$\frac{\partial }{\partial t}\left({\rho }_g{\alpha }_g+{\rho }_l{\alpha }_l\right)+\frac{\partial }{\partial x}\left({\rho }_g{{\alpha }_{g}v}_g+{\rho }_l{{\alpha }_{l}v}_l\right)=0$$

(1)

ρ_g andρ_l are the densities of the gas phase and liquid phase, respectively. v_g and v_l are the velocities of the gas phase and liquid phase, respectively. α_g and α_l are the cross-sectional gas content and cross-sectional liquid content.

Momentum conservation equation

$$\frac{\partial }{\partial t}\left({\rho }_g{{\alpha }_{g}v}_g+{\rho }_l{{\alpha }_{l}v}_l\right)+\frac{\partial }{\partial x}\left(P+{\rho }_g{\alpha }_g{v}_g^2+{\rho }_l{{\alpha }_{l}v}_l^2\right)=-{\rho }_m{v}_{m}g\sin \theta -F_f$$

(2)

P is the wellbore pressure, ρ_m and v_m are the density and velocity of the gas-liquid mixture, respectively, f is the wellbore friction coefficient, and D is the wellbore diameter.

Energy conservation equation

$$\frac{\partial }{\partial t}\left[{\rho }_g{\alpha }_g\left(e_g+\frac{1}{2}{v}_g^2\right)+{\rho }_l{\alpha }_l\left(e_l+\frac{1}{2}{v}_l^2\right)\right]+\frac{\partial }{\partial x}\left({\rho }_g{{\alpha }_{g}v}_g{h}_{Tg}+{\rho }_l{{\alpha }_{l}v}_l{h}_{Tl}\right)+\frac{Q}{\pi r_i^2}=0$$

(3)

e_g and e_l are the specific internal energies of the gas and liquid phases, respectively. h_Tg and h_Tl are the total specific enthalpies, defined as the sum of specific enthalpy, kinetic energy, and potential energy. Q represents the energy change due to heat exchange between the formation and the wellbore. r_i is the wellbore radius.

$$e_g=h_g-\frac{P}{{\rho }_g}$$

(4)

$$e_l=h_l-\frac{P}{{\rho }_l}$$

(5)

$$h_{Tg}=h_g+\frac{v_g^2}{2}-gx$$

(6)

$$h_{Tl}=h_l+\frac{v_l^2}{2}-gx$$

(7)

Where h_g and h_l represent the specific enthalpies of the gas phase and liquid phase, respectively.

2.2. Constitutive equations

Due to the large number of unknowns in these conservation equations, it is necessary to supplement them with constitutive equations. In geothermal wellbore modeling, considerations include the slip between gas and liquid phases, friction with the wellbore wall, and heat transfer between the wellbore and the surrounding formation.

2.2.1. Mixture equation

The relationship between the cross-sectional gas content and the mass liquid content is as follows:

$${\alpha }_l+{\alpha }_g=1$$

(8)

The mixture density is described as:

$${\rho }_m={\alpha }_g{\rho }_g+{\alpha }_l{\rho }_l$$

(9)

2.2.2. Wellbore friction

In geothermal drilling, the gas-liquid mixture flowing upward along the wellbore causes a significant frictional pressure drop. The expression for F_f is as follows:

$$F_f=\frac{f{\rho }_m{v}_m\left|v_m\right|}{2D}$$

(10)

The f in Equation (7) represents the wellbore friction coefficient, which can be expressed as:

$$f=\left\{\begin{array}{cc}\frac{16}{Re}\hfill & Re\le 2400\hfill \\ \frac{1}{16}{\left\{log\left[\frac{\xi }{3.7r_i}-\frac{5.02}{Re}log\left(\frac{\xi }{3.7r_i}+\frac{13}{Re}\right)\right]\right\}}^{-2}\hfill & Re>2400\hfill \end{array}\right.$$

(11)

Where Re is the Reynolds number, ξ is the wellbore roughness. In this paper, the roughness is assumed to be 0.00175 millimeters. The Reynolds number is calculated as follows:

$$Re=\frac{{\rho }_m\left|v_m\right|2r_i}{{\mu }_m}$$

(12)

In Equation (9), μ_m represents the viscosity of the gas-liquid mixture, and its expression is as follows:

$${\mu }_m={\alpha }_g{\mu }_g+{\alpha }_l{\mu }_l$$

(13)

2.2.3. Formation-Wellbore Heat Exchange

In geothermal drilling, there is heat transfer between the formation and the wellbore. Heat is transferred from the formation through the cement and casing into the wellbore, and the fluid inside the wellbore is transferred back to the formation through the same path. Here, we use a heat transfer equation developed by Hasan [14], applicable to any wellbore, which is expressed as follows:

$$Q=\frac{2\pi r_{ti}{U}_t{k}_f\left(T_{ei}-T_f\right)}{k_f+\left(r_{ti}{U}_t{T}_D\right)}$$

(14)

k_f is the formation thermal conductivity, T_ei is the formation temperature unaffected by wellbore heat transfer, and T_f is the fluid temperature inside the drill string. Radial heat exchange occurs between the formation and the wellbore fluid. When heat is transferred from the formation to the wellbore, it first passes through the cement, then the casing, and finally through the outer wall of the drill string into the fluid inside the drill string. Each of these parts has its own thermal resistance, so the overall heat transfer coefficient U_t can be considered as the series combination of these individual resistances, as represented by Equation (16), T_D is a temperature function related to the dimensionless time t_d.

$$\frac{1}{U_t}=\frac{1}{h_{ti}}+\frac{r_{ti}ln\left(\frac{r_{to}}{r_{ti}}\right)}{k_t}+\frac{r_{ti}}{r_{to}\left(h_c+h_r\right)}+\frac{r_{ti}ln\left(\frac{r_{co}}{r_{ci}}\right)}{k_c}+\frac{r_{ti}ln\left(\frac{r_{wb}}{r_{co}}\right)}{k_{cem}}$$

(15)

$$T_D=ln\left[e^{-0.2t_d}+\left(1.5-0.3719e^{-t_d}\right){t_d}^{\frac{1}{2}}\right]$$

(16)

2.2.4. Slip Model

In the wellbore, the gas moves faster than the liquid, often resulting in slip. Since Zuber and Findlay introduced the drift-flux model, many researchers have studied various drift-flux models [15,16,17,18,19]. Investigations have shown that among these models, the one developed by Shi et al [10] has better applicability. Therefore, this paper adopts the drift-flux model established by Shi et al [10]. The two-phase velocity formula is described as：

$$v_g={C_{0}v}_m+v_d$$

(17)

Where v_m is the mixture velocity, and v_d is the drift velocity.

$$v_m=v_{sg}+v_{sl}={\alpha }_g{v}_g+{\alpha }_l{v}_l$$

(18)

$$v_d=\frac{m\left(\theta \right)\left(1-C_0{\alpha }_g\right)K\left({\alpha }_g\right)v_c}{C_0{\alpha }_g\sqrt{\frac{{\rho }_g}{{\rho }_l}}+1-C_0{\alpha }_g}$$

(19)

In this slip model, m(θ) is a function of the wellbore inclination angle, C₀​ is the profile coefficient, K is a function of the cross-sectional gas content, and v_c is the characteristic velocity.

In a gas-water dominated system, the expression for m(θ) is:

$$m\left(\theta \right)=1.85{\left(cos\left(\theta \right)\right)}^{0.21}{\left(1+sin\left(\theta \right)\right)}^{0.95}$$

(20)

Here, the improved K by Pan et al. [17]. based on Shi's model [10] is used.

$$K\left({\alpha }_g\right)=\left\{\begin{array}{c} 1.53 {\alpha }_g\le {\alpha }_1 \\ 1.53+0.5\left(C_0{K}_u-1.53\right)\left(1-cos\left(\pi \frac{{\alpha }_g-{\alpha }_1}{{\alpha }_2-{\alpha }_1}\right)\right) {\alpha }_1<{\alpha }_g< {\alpha }_2 \\ C_0{K}_u { \alpha }_g\ge {\alpha }_2\end{array}\right.$$

(21)

where K_u is the "critical Kutateladze number" related to the pipe diameter.

$$K\left(\widehat{D}\right)={\left(\frac{g\left({\rho }_l-{\rho }_g\right)}{{\sigma }_{gl}}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}D$$

(22)

The expression for the characteristic velocity v_c is:

$$v_c={\left(\frac{{\sigma }_{gl}\left({\rho }_l-{\rho }_g\right)}{{{\rho }_l}^2}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$$

(23)

The formula for the profile parameter C₀ is:

$$C_0=\frac{A}{1+\left(A-1\right){\gamma }^2}$$

(24)

This equation ensures that C₀ has a constant value during bubble flow or slug flow, and equals 1 during annular mist flow. The expression for the parameter γ is:

$$\gamma =\frac{\beta -B}{1-B}$$

(25)

$$\beta =max\left({\alpha }_g,F_v\frac{{\alpha }_g\left|v_m\right|}{v_{sgf}}\right)$$

(26)

F_v is a constant used to adjust the sensitivity of C₀ to velocity. v_sgf is the critical velocity value that prevents the liquid film from falling back with the gas flow.

$$v_{sgf}=K_u{\left(\frac{{\rho }_l}{{\rho }_g}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{v}_c$$

(27)

After experiments conducted by Shi et al., [10] it was found that in the gas-water system, the optimal parameter A is equal to 0, α₁ = 0.06, and α₂ = 0.21. Therefore, the calculated value of C₀ is equal to 1.

3. Numerical methods

Finite Difference Method (FDM) is a common numerical analysis technique used to solve differential equations. It discretizes the differential equations in a specified region, transforming continuous problems into discrete ones. It has been widely applied in wellbore modeling [20,21]. Due to the highly nonlinear and strongly coupled nature of the geothermal wellbore model, analytical solutions are not feasible. Therefore, numerical methods are needed to approximate the desired solutions. The flowchart below illustrates the solution process in this paper.

• Initialize Pressure and Temperature: First, set the initial pressure and temperature values, then divide the wellbore into grids.
• Calculate Velocity and Phase Fractions: Use the mass conservation equation and momentum conservation equation to determine the fluid velocities and phase fractions. This step involves calculating the movement properties of the fluids within the wellbore.
• Calculate Temperature: Apply the energy conservation equation to compute the temperature, considering the heat conduction within the wellbore fluids.
• Newton-Raphson Iteration: Since the discretized equations form a nonlinear system, use the Newton-Raphson iteration method. Based on the current values of velocity, phase fractions, temperature, and pressure, compute the Jacobian matrix and iterate to solve the nonlinear equations.
• Error and Convergence Check: Compare the calculated residuals with the convergence criteria. If the error is less than the threshold of 10^-6, the solution is considered converged, and the loop can end.
• Update Solution and Continue Loop: If the solution has not converged, update the velocity, phase fractions, temperature, and pressure values based on the Newton-Raphson iteration results and error assessment. Then, return to step 1 to continue the iteration process.
• Advance to Next Time Step: Once the simulation converges, proceed to the next time step and continue simulating the wellbore conditions at the next moment.

3.1. Discretization of Equations

3.1.1. Grid Division

To accurately simulate the geothermal wellbore model, the wellbore is divided into i_n discrete grid cells along its length, with a total length of L and an increment of Δx for each segment (see Figure 2). Detailed calculations of pressure, temperature, velocity, and phase fractions are performed from the initial grid point at the wellhead, iterating down to the bottom of the well.

3.1.2. Mass conservation equation

The mass conservation equation is discretized in time and space using a forward difference scheme, given as follows:

$$\frac{{\left({\rho }_g{\alpha }_g+{\rho }_l{\alpha }_l\right)}_i^{n+1}-{\left({\rho }_g{\alpha }_g+{\rho }_l{\alpha }_l\right)}_i^n}{∆t}+\frac{{\left({\rho }_g{{\alpha }_{g}v}_g+{\rho }_l{{\alpha }_{l}v}_l\right)}_{i+1}^n-{\left({\rho }_g{{\alpha }_{g}v}_g+{\rho }_l{{\alpha }_{l}v}_l\right)}_i^n}{∆x}=0$$

(28)

3.1.3. Momentum conservation equation

The momentum conservation equation is discretized as follows:

$$\begin{gathered} \frac{{\left({\rho }_g{{\alpha }_{g}v}_g+{\rho }_l{{\alpha }_{l}v}_l\right)}_i^{n+1}-{\left({\rho }_g{{\alpha }_{g}v}_g+{\rho }_l{{\alpha }_{l}v}_l\right)}_i^n}{∆t}+\frac{{\left(P+{\rho }_g{\alpha }_g{v}_g^2+{\rho }_l{{\alpha }_{l}v}_l^2\right)}_{i+1}^n-{\left(P+{\rho }_g{\alpha }_g{v}_g^2+{\rho }_l{{\alpha }_{l}v}_l^2\right)}_i^n}{∆x} \\ =-\left(\frac{{\left({\rho }_m{v}_m\right)}_{i+1}^n-{\left({\rho }_m{v}_m\right)}_i^n}{∆x}\right)g\sin \theta -\frac{f}{2D}\left(\frac{{\left({\rho }_m{v}_m\left|v_m\right|\right)}_{i+1}^n-{\left({\rho }_m{v}_m\left|v_m\right|\right)}_i^n}{∆x}\right) \end{gathered}$$

(29)

3.1.4. Energy conservation equation

The Energy conservation equation is discretized as follows:

$$\begin{gathered} \frac{{\left[{\rho }_g{\alpha }_g\left(e_g+\frac{1}{2}{v}_g^2\right)+{\rho }_l{\alpha }_l\left(e_l+\frac{1}{2}{v}_l^2\right)\right]}_i^{n+1}-{\left[{\rho }_g{\alpha }_g\left(e_g+\frac{1}{2}{v}_g^2\right)+{\rho }_l{\alpha }_l\left(e_l+\frac{1}{2}{v}_l^2\right)\right]}_i^n}{∆t} \\ +\frac{{\left({\rho }_g{{\alpha }_{g}v}_g{h}_{Tg}+{\rho }_l{{\alpha }_{l}v}_l{h}_{Tl}\right)}_{i+1}^n-{\left({\rho }_g{{\alpha }_{g}v}_g{h}_{Tg}+{\rho }_l{{\alpha }_{l}v}_l{h}_{Tl}\right)}_i^n}{∆x}+\frac{Q}{\pi r_i^2}=0 \end{gathered}$$

(30)

3.2. Boundary conditions

Boundary conditions in the geothermal wellbore model are crucial for the model's accuracy and need to be precisely defined.

3.2.1. At the wellhead

The temperature at the wellhead is determined by the temperature of the drilling fluid at the wellhead, and is thus expressed as:

$$T_h\left(0,t\right)=T_m$$

(31)

The pressure at the wellhead is the sum of the dynamic pressure and the backpressure, thus the expression for the pressure at the wellhead is:

$$P\left(0,t\right)=P_C+\frac{1}{2}{\rho }_m\left(t,x\right){\left(v_m\left(t,x\right)\right)}^2$$

(32)

3.2.1. At the bottom-hole

Due to the significant influence of the formation temperature at the bottom of the wellbore, the temperature inside the wellbore is expressed as Equation (34), where ΔT represents the temperature loss, determined according to Fourier's law of heat conduction, and ΔZ represents the distance between the wellbore and the formation.

$$T\left(L,t\right)=T_{ei}-∆T$$

(33)

$$∆T=\frac{Q∆Z}{k_f}$$

(34)

Assume the bottom-hole pressure is equal to the formation pressure.

$$P\left(L,t\right)=P_f$$

(35)

3.3. Newton-Raphson Iteration

After discretization, the conservation equations remain nonlinear. Therefore, the Newton-Raphson method is used for iterative solution. The Newton-Raphson method has many advantages, such as efficiency and fast convergence, and it has shown good application effects in wellbore modeling [22,23,24].

Assuming the solution vector is:

$$\mathit{u}=u{\left({\rho }_g,{\alpha }_g,v_g,P,T\right)}^T$$

(36)

The residual vector F is calculated from the conservation equations. The mass conservation equation, momentum conservation equation, and energy conservation equation are represented as F₁, F₂, and F₃ respectively. Additionally, the equations from the drift-flux model and the heat transfer model are used as auxiliary residual equations F₄ and F₅.

$$\mathit{F}\left(\mathit{u}\right)=\left(\begin{array}{c}{F}_1\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_2\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_3\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_4\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_5\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\end{array}\right)$$

(37)

The increment Δu is obtained by multiplying the inverse of the Jacobian matrix by the residual vector.

$$∆\mathit{u}={\left[\begin{array}{ccccc}\frac{\partial F_1}{\partial {\rho }_g}& \frac{\partial F_1}{\partial {\alpha }_g}& \frac{\partial F_1}{\partial v_g}& \frac{\partial F_1}{\partial P}& \frac{\partial F_1}{\partial T}\\ \frac{\partial F_2}{\partial {\rho }_g}& \frac{\partial F_2}{\partial {\alpha }_g}& \frac{\partial F_2}{\partial v_g}& \frac{\partial F_2}{\partial P}& \frac{\partial F_2}{\partial T}\\ \frac{\partial F_3}{\partial {\rho }_g}& \frac{\partial F_3}{\partial {\alpha }_g}& \frac{\partial F_3}{\partial v_g}& \frac{\partial F_3}{\partial P}& \frac{\partial F_3}{\partial T}\\ \frac{\partial F_4}{\partial {\rho }_g}& \frac{\partial F_4}{\partial {\alpha }_g}& \frac{\partial F_4}{\partial v_g}& \frac{\partial F_4}{\partial P}& \frac{\partial F_4}{\partial T}\\ \frac{\partial F_5}{\partial {\rho }_g}& \frac{\partial F_5}{\partial {\alpha }_g}& \frac{\partial F_5}{\partial v_g}& \frac{\partial F_5}{\partial P}& \frac{\partial F_5}{\partial T}\end{array}\right]}^{-1}\left(\begin{array}{c}{F}_1\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_2\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_3\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_4\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\\ F_5\left({\rho }_g,{\alpha }_g,v_g,P,T\right)\end{array}\right)$$

(38)

The iterative formula for solving the solution vector in conjunction with the auxiliary equations mentioned earlier is given by Equation (40).

$${\mathit{u}}_{\mathit{i}+1}={\mathit{u}}_{\mathit{i}}+∆\mathit{u}$$

(39)

4. Model verification and analysis

To validate the effectiveness of the model, three sets of comparative experiments were conducted using geothermal wells from different locations. The first set involved a geothermal well in the Rongxi area of Xiongan, China. The parameters of this well are listed in Table 1. The second set of experiments used the SNLG87-29 well data from Garg et al. [25]. Finally, the third set utilized data from the No. 6 well at point D in northwest Iran, as mentioned in the article by Zolfagarroshan et al. [26]. First, the parameters of the well in the Rongxi area of Xiongan are listed in Table 1.

Based on the data in the table, calculations using the model produced results as shown in Figure 3. Figure 3-a presents the measured temperature data inside the drill pipe. After comparing and validating with the models by Akbar et al. and Tonkin et al., the results are as follows: all three wellbore models provide accurate temperature calculations. The mean relative error (MRE) for this model is 1.323%. For Akbar's model, the MRE is 1.999%, and for Tonkin et al.'s model, the errors are 1.78%. This indicates that the temperature calculation model in this study is relatively accurate.

Figure 3-b shows the trend of bottom hole pressure over time. As the drilling depth increases over time, the bottom hole pressure increases. However, due to various uncertainties in the formation, the bottom hole pressure may suddenly change, necessitating adjustments to the surface backpressure to control the bottom hole pressure. The figure shows that the model in this study provides a good fit during the stable drilling phase, outperforming the other two models. When the bottom hole pressure suddenly changes, the MRE in the calculated results is 1.491%, which is smaller than the MRE calculated by the other two models.

To clearly and intuitively demonstrate the pressure calculation results of various geothermal wellbore models, pressure-depth curves were plotted at different time points: t=0, 400, and 800.

Figure 5 shows that both the wellhead pressure and the bottom hole pressure increase with time. The pressure-depth curves calculated by this model at different time points match well with the actual measured data. The Mean Relative Errors (MRE) at the three different time points are 1.254%, 1.254%, and 1.339, respectively.

Figure 6 shows the pressure results calculated by Akbar's model and Tonkin's model. The Mean Relative Errors (MRE) calculated by Akbar's model at different time points are 2.036%, 2.035%, and 2.025%, while the MREs calculated by Tonkin's model at different time points are 2.037%, 2.029%, and 2.029%.

The SNLG87-29 well was designed to match the temperature and pressure data recorded in 1993. The parameters are listed in Table 2.

Based on the existing parameters and some initial parameter assumptions, calculations were performed for this geothermal well. The results are shown in Figure 6.

Figure 6-a illustrates the pressure variation with depth. It can be clearly observed that the measured data align well with the model presented in this study, with an MRE of: 1.824%. Given the shallow depth and low pressure of the well, all three models show good fitting results. The MRE for Akbar's model is 4.624%, and for Tonkin's model, it is: 3.213%.

Figure 6-b depicts the temperature variation with depth. In the initial 75 meters, all three models fit the measured data well. However, as the temperature gradient changes, the temperature variation with depth decreases. Akbar's model retains the trend of the original data curve but has a higher MRE, which is: 2.203%. The model presented in this study shows a similar trend to Tonkin's model, but with an MRE of: 1.209%, which is lower than Tonkin's.

This indicates that the model presented in this study performs well even in high-temperature, low-pressure wells.

Zolfagharroshan et al. [26] mentioned detailed parameters of Well No. 6, listed in Table 3.

Based on the parameters of Well No. 6 listed in Table 3, a comparison of the results calculated by the three models with the measured data is shown below:

Figure 7(a) shows that the pressure-depth curve of this well has two inflection points, where the rate of pressure change with depth significantly varies. The model proposed in this paper closely matches the measured data, with a calculated Mean Relative Error (MRE) of 2.937%. The Akbar model and the Tonkin model also exhibit good pressure curves but show larger errors compared to the measured data and lag in the pressure increasing trend. The calculated MREs are 6.55% and 3.185%, respectively.

Figure 7 illustrates that the temperature calculation results of the three models have relatively small errors. However, this model has an MRE of 0.912%, which is lower than the other two models. Meanwhile, the temperature calculation results of the Akbar model are superior to those of the Tonkin model, with MREs of 1.553% and 1.814%, respectively.

5. Conclusions

This paper presents a geothermal wellbore model that demonstrates strong adaptability under three different temperature and pressure well conditions. The model integrates conservation equations with the drift-flux model and heat transfer model. After discretizing the equations, the Newton-Raphson iteration method is used to calculate pressure and temperature. Validation results show that the model performs well in three different wellbore environments, achieving high accuracy by taking into account multiple factors affecting heat transfer. Experimental comparison results indicate that the model has a low Mean Relative Error (MRE), demonstrating its high precision and reliability in calculating bottom hole pressure and temperature. Overall, this model provides reliable and accurate results, making it a valuable tool for optimizing drilling parameters and ensuring drilling safety in the geothermal energy field.