Modeling of two-phase flow of high temperature geothermal production wells in the Yangbajing geothermal field, Tibet

Two-phase flow (flow of water in both liquid and gas phase, containing a non-condensable gas such as CO2) in the wellbore is one of the most important processes in the production performance and wellbore scaling evaluations of high temperature geothermal wells. This paper first describes the discharge tests in the Yangbajing geothermal field, Tibet. Next, a simple model for governing the two-phase flow in the presence of CO2 in the wellbore is constructed and a robust calculation method is presented. The model is applied to three production wells in the Yangbajing geothermal field. The results show that the velocity difference between the gas and liquid phase should be included in the model. Ignoring this difference (i.e., homogeneous model) would result in a significant deviation between measurements and calculations. A drift flux model (DFM) describes the velocity difference, where the specifics of the particular model can have significant effects on the results for the pressure and temperature profiles in the wellbore. Three commonly used DFMs are compared to estimate their performance. The calculated wellhead pressure and temperature are in the range of 2,3 bar and 125°C–135°C for all three wells at a production rate of about 20 kg/s. The estimated wellhead gas mass fraction is between 3% and 8%. Considering CO2 content, three different scenarios were evaluated, although the effect on the pressure and temperature profiles were limited. However, CO2 content has a much more significant effect on the flash depth, which is an important parameter for the estimation of calcite scaling.


Introduction
As a natural resource with low carbon content and huge reserves, geothermal energy has become a promising renewable energy source and has started to attracted more attention. Theoretical calculations show that the energy reserves in the upper 10 km of the Earth's crust are approximately 3.6 ✕ 10 14 GWh (Lund et al., 2008), about 2.17 million times the global energy consumption of 2012. There are roughly 88 countries (regions) making use of geothermal energy directly (Lund and Toth, 2021), and more than 30 countries (regions) using geothermal energy for power generation (Huttrer, 2021). Most of the geothermal energy is used by geothermal heat pumps. About 12.9% is used for electricity generation (Lund et al., 2022), accounting for just 0.3% of the total power generation and 1.5% of the power generated from renewables (Zhu et al., 2015). The United States is the leading country in geothermal applications. It has the largest installed capacity of geothermal power production, reaching 28.8% of the total geothermal power production. Additionally, it has the second largest installed capacity of direct geothermal utilization. China is number one when it comes to direct geothermal utilization, reaching 25.2% of the total. However, the installed capacity of geothermal power generation is only 27.28 MWe, accounting for just 0.2% of the total geothermal power production and ranking 18th in the world (China Geological Survey of Natural Resources Ministry et al., 2018).
China has rich geothermal resources, with a value of 3.06 ✕ 10 18 kWh/a accounting for 7.9% of the total global geothermal energy reserves (Jiang, 2012). However, the distribution of geothermal resources is uneven but with obvious regularity and regionality. At present, a basic primary pattern for geothermal energy utilization in China has been formed (Zhu et al., 2015). Yangbajing in Tibet features geothermal power generation. Geothermal heating is widely used in Tianjin, Xi'an and Beijing. Groundwater and sea water source heat pumps for heating and cooling are widely used in Chongqing and Dalian. The coastal regions in South and East China are popular for recuperation and tourism related utilization such as spa hotels and hot spring resorts. The first high-temperature geothermal power generation station in China was built in Yangbajing, Tibet in 1977. The total installed capacity is currently 26.18 MW. In the past few decades, the power station has contributed 50% of Lhasa's summer power supply and 60% of winter power supply (Zhang et al., 2019). However, the operation of the power station was suspended due to the low power generation efficiency, poor economic efficiency (Zhang et al., 2019), geothermal wastewater (Guo and Wang, 2009) and scaling problems (Zhou, 2003).
The wellbore is the only conduit available to extract geothermal fluid from reservoirs. Some important physical and geochemical processes such as phase change, gas exsolution and mineral precipitation are likely to occur as a result of the fluid flowing upward, and these processes can significantly affect the production performance of a wellbore. Two-phase flow with phase change is generally the basis for another processes in wellbores. Due to the difficulty in measuring directly some important parameters as the harsh operating conditions in geothermal wellbores prevent the installation of measuring equipment. Therefore, wellbore flow simulator is an invaluable tool for obtaining the information on flow behavior. Barelli et al. (1982) constructed a model that simulates a twophase flow in geothermal wells that contain CO 2, and gives an estimate of temperature and pressure profiles which are then compared to measurements of the specific well. Bjornsson et al. (1993) developed the geothermal wellbore simulator HOLA. The simulator can reproduce the measured flowing temperature and pressure profiles in flowing wells and determines the relative contribution of each feed zone for a given discharge condition. Pan and Oldenburg. (2014) originally developed the simulator T2Well in the TOUGH2 framework to model two-phase flows in coupled wellbore-reservoir systems. Vasini et al. (2018) plugged EWASG (Battistelli et al., 1997) into T2Well in order to extend the functionality of computing two-phase flow for high enthalpy geothermal systems. The module EWASG can consider fluids containing dissolved solids and one non-condensable gas (NCG) such as CO 2 , CH 4 , H 2 S, H 2 or N 2 . Gunn and Freeston. (1991) developed an integrated geothermal wellbore simulator called WELLSIM, which can perform calculations including discharge tests, fluid composition and properties, deliverability curve prediction, statistical and graphical matching analysis of downhole pressure and temperatures and downhole measurement analysis. Hasan and Kabir. (2010) presented a robust model for a two-phase flow in geothermal wells using the drift-flux approach. Deendarlianto and Itoi. (2021) constructed a wellbore model to numerically investigate the effects of CO 2 gas on the two-phase flow characteristics in geothermal production wells. Lei et al. (2022) evaluated the two-phase flow in the geothermal wells by the Shi et al. (2005) model. Li et al. (2020) used HOLA and WELLSIM to simulate the two-phase flow of a high temperature geothermal well in the Kangding geothermal field for quantitative assessment of calcite scaling. Recently, Tonkin et al. (2021) reviewed geothermal wellbore models and compared their differences and effects on simulations. Furthermore, Canbaz et al. (2022) performed a comprehensive literature review of studies about wellbore dynamics involving geothermal fluids' physical and thermodynamic behavior during production and injection in the presence of CO 2 along the wellbore.
This paper aims to estimate the two-phase flow and related production performance in the Yangbajing geothermal field by comparing numerical simulation with those obtained from discharge testing. It is organized as follows: Section 2 briefly introduces the field condition, and pressure and temperature measurements during discharge testing. Section 3 describes the details of the wellbore model. In Sections 4 and 5 the proposed wellbore model is used to simulate the two-phase flow of typical production wells in the Yangbajing geothermal field. Finally, the main observations and conclusions are summarized in Section 6.
2 Study area

Geological settings
The Yangbajing geothermal field is located in a fault basin of the Nyainqentanglha range, about 90 km northwest of Lhasa, at an elevation of 4,290-4,500 m with a high altitude in the northwest and a low altitude in the southeast. There are strong tectonic activities in the geothermal field, with three groups of faults developing in northeast, northwest and north-south directions. The China-Nepal highway divides the field into two parts, the southern and northern. The south area is composed of porous Quaternary strata, while the north area is composed of a combination of porous Quaternary strata and Himalayan fractured granite (Duo, 2003). The thermal groundwater system consists of two reservoirs at different depths. The shallow geothermal reservoir is located 180-280 m below the surface, and the lithology is composed of Quaternary alluvial sandy gravel layers, moraine gravel and weathered granite crust on top of the bedrock. The cap rock of the geothermal reservoir is composed of muddy gravel or silty clay layers with varying thickness. The bedrock at the bottom of the geothermal reservoir is composed of early Himalayan granite and tuff, and mylonite granite is found locally in the northern part of the geothermal field (Duo, 2003;Guo et al., 2007).

Frontiers in Earth Science
frontiersin.org 02 According to the magnetotelluric exploration results, there is a low resistance layer with a resistivity of 5 Ω m at a depth of 5 km in the northern part of the Yangbajing geothermal field, which is inferred to be a cooling magma chamber and the heat source of the geothermal field (Chen, et al., 1996;Duo, 2003). A fault fracture zone is present in the piedmont of the north, with a local width of about 1 km. The shallow normal faults are connected with the deep strike-slip faults, effectively connecting the deep heat source and the shallow reservoir. The results

FIGURE 1
The location and geological profile of the study area:(A) plan, (B)cross section (Modified from Duo, 2003;Xu et al., 2018).

Frontiers in Earth Science
frontiersin.org 03 of analysis of the hydrogen and oxygen isotopes of the groundwater in the thermal field (Guo et al., 2010) indicates that the subterranean hot water in this area is the result of a large amount of melt water from the Nyainqentanglha Mountain flowing into the ground along the fault zone, circulating deep, and being heated by the deep magma heat source (Duo, 2003) (Figure 1).
The fluid is mainly in liquid state containing a little gas. Guo. (2012) reviewed the hydrogeochemistry of the Yangbajing geothermal field, indicating that 1) the TDS of the hot fluid is about 1,266 mg/L and Na + , Cl − and HCO 3 -are the dominant ions, and 2) CO 2 originated from the metamorphism of the Nyenchen Tonglha core complex is a major component of geothermal gases. Gas sampling by a separator at the wellhead was carried out by Zhao et al. (1998). The results showed that CO 2 content in the gas-phase was between 21.4 and 66.10 mmol/kg. CO 2 was the dominant composition of non-condensable gas ( Table 1).
The geological survey in the Yangbajing geothermal field began in 1976. During the exploration stage, 45 exploration wells were drilled and a thermal field with an area of 42 km 2 was delineated. In the 1980s, another 25 production wells and 13 recharge wells were successively drilled. Since September 1977, the first 1 MW test unit has successfully been generating power. The power generation capacity had been expanded to 25.18 MW by 1991. By the end of 2017, the field consisted of four production wells in the south, 12 production wells in the north, and five recharge wells, reaching a total installed capacity of 24.18 MW and a cumulative power generation of more than 33×10 8 kW h. The fluid production rate was about 12,000 t/d in the peak production period. The initial fluid temperature of most production wells was 125°C-145°C, with a pressure of 1.76-3.72 bar. The temperature in the shallow geothermal reservoir is ranging from 140°C to 160°C. The pressure in the reservoir varies from about 0.8 MPa to 3.0 MPa, depending on the depth. However, the pressure and temperature of the production wells showed a decreasing trend with time before 2012. Pressure and temperature were gradually restored after recharge normalization after 2012 (Xu et al., 2018). The long-term fluid extraction leaded to considerable surface subsidence around the geothermal fields (Li et al., 2015).

Discharge testing
In order to analyze the flow of the wellbore, the discharge testing was carried out from April to May 2020 in three typical wells, namely, ZK303, ZK324, and ZK304. The ZK303 well has a depth of 300 m, where the first 142.5 m in depth is in the Quaternary layer. Below 142.5m, the well extends into the fractured layer. The maximum depth of mechanical descaling of this well is about 150 m, and the scaling occurs above~80 m ( Figure 2A). The ZK324 well has a depth of 90.13m, and the entire well section is in the Quaternary sandy gravel layer. The maximum depth of mechanical descaling of this well is about 86 m, and the scaling occurs above~70 m ( Figure 2B). The ZK304 well has a depth of 206.54 m, with the Quaternary porous layer from 0 m to 142 m in depth and the fractured layer below the depth 142 m. The maximum depth of mechanical descaling is 117.17 m, and the scaling occurs between 80 m and 130 m ( Figure 2C).
The James lip pressure method was adopted in the discharge testing, using the Kuster K10 PT Geothermal Instrument as the main test tool. Discharge stimulation was conducted after the static PT test (i.e., without discharge) and a continuous PT test was performed from the bottom to the top with wellhead fluid discharge. The flow rate was measured at the wellhead. Figure 3 shows the test results. It can be seen that the temperature and pressure at the bottom of the well are not much different under the two testing conditions, but the difference at the upper part of the wells is relatively large. In the static state, the water level is below the wellhead elevation; the water levels of ZK303 and ZK304 are at the depth of 50-60 m, and the water level of ZK324 is at the depth of about 10 m. Due to the combined effects of atmospheric environment, geothermal reservoirs and surrounding rocks near the wellbore, the temperature shows a large gradient. In the dynamic state, the wellhead temperature is relatively low due to the influence of atmospheric environment, but the temperature at the depth of 10 m is 118°C-135°C.
To further analyze the phase characteristics from the downhole to the wellhead under dynamic condition, the data of temperature and pressure are plotted on the water phase diagram. It can be seen from Figure 4 that the fluid in the three wells is in a single liquid-phase state at the bottom. As the fluid flows upward, the pressure gradually decreases and the fluid begins to evaporate when it drops to the saturation pressure of the corresponding temperature. As a result, the temperature decreases and the system goes to two-phase state. In the two-phase state, the relationship between temperature and pressure is locked by the water saturation line. Due to the influence of the surface environment, ZK324 and ZK304 experienced strong condensation at the wellhead during the testing, which caused the fluid state to return to single liquidphase.
The phase diagram ( Figure 4) shows that the change in pressure and temperature is almost along the path of the saturation line of pure water. It indicates that CO 2 in the system is limited. The measurements of wellhead steam compositions (Zhao et al., 1998) showed that the CO 2 mass weight fractions in steams are 0.29% and 0.12% for ZK303 and ZK304, respectively. Furthermore, a large amount of calcite scaling was found in the wellbore during several years of operation (Zhou, 2003), which is strongly associated with the presence of CO 2 . Therefore, the presence of CO 2 should be considered in the evaluation of twophase flow.
3 Simulation method 3.1 Flow with phase change in onedimension wellbores The components of fluid in the Yangbajing geothermal fields contain H 2 O and CO 2 and phase change may occur when the fluid flows upward in the wellbore. Based on mass, momentum and energy conservations, a one-dimensional non-isothermal flow model with phase change can be constructed for the wellbore. Assuming steady state and uniform distribution of variables (i.e., temperature, pressure, gas-phase fraction, density and velocity) along the wellbore cross-section and ignoring the axial heat conduction of the fluid in the wellbore, the equations can be described by the following form (Tonkin et al., 2021) Total mass: where.
-P is the total pressure comprised of P H2O + P CO2 .
-α is the void fraction defined as the fraction of the channel cross sectional area that is occupied by gas phase.
-ρ l and ρ g are the liquid and gas densities, respectively. -u l and u g are the liquid and gas velocities, respectively. -x CO2 l and y CO2 g are the CO 2 mass fractions in liquid and gas phase, respectively. -_ m l and _ m g are the mass rates of liquid and gas phase, respectively. _ m l and _ m g are defined as respectively (1 − x) _ m and x _ m, where _ m is the total mass rate. -h l and h g are the enthalpies of liquid and gas phases, respectively. -A is the wellbore cross-section area.

Frontiers in Earth Science
frontiersin.org 05 -z is the elevation -d is the diameter -g is the gravitational acceleration.
-f is the friction factor as a function of Reynolds number Re. When Re Re )] μ m and ε are the mixture viscosity and roughness of the wellbore, respectively. -ρ m is the density of mixture as ρ m (1 − α)ρ l + αρ g.
-u m is the velocity of mixture as . -q is a term that describes the surrounding rock heat exchange, defined as q 2πr 1 U(T f − T wb ). With the assumption that only heat conduction in the infinite surrounding rock is considered, the parameter U u r1uf(tD)+kT depends on wellbore radius (r 1 ), completion radius (r 2 ), thermal conductivity of the formation (k T ) and the wellbore completion (k cem ), density (ρ r ) and specific heat (c r ) of the formation, and time (t), where u kcem r1 ln ( r 2 r 1 ) ,

FIGURE 3
Temperature and pressure test results, (A) pressure, (B) temperature.

FIGURE 4
Phase behavior along wellbore under dynamic test condition.
With the assumption that momentum related terms in Eqs 3, 4 are negligible, and by considering the boundary conditions at the downhole, Eqs 1-4 can be reduced to The mass rate _ m and the CO 2 mass fraction γ are often known from the field measurements. For a single-phase flow, only Eqs 5 and 7, 8 are needed, as the CO 2 mass fraction in liquid phase x CO2 l is equal to a constant value of γ. The three variables of P, h and u are used as the primary variables and all the other parameters are functions of these three variables. For a two-phase flow, the unknown variables are ρ l , ρ g , u l , u g , α, x CO2 l , y CO2 g , P, h l , h g and gas mass fraction x. Therefore, additional equations are required. When the fluid is in a two-phase state, the partial pressure of water is dependent on temperature as Although the uneven distribution of gas bubble would result in the non-equilibrium of CO 2 between the gas and liquid phase, this paper also assumes that the component of CO 2 has enough chance to reach an equilibrium state between the two phases. According to Henry's law, the relation between the partial pressure of CO 2 and CO 2 mass fraction in liquid phase can be expressed as where H(T) is Henry's law coefficient, only considering depending on temperature. This assumption is reasonable for low CO 2 partial pressure (less than 10 bar). For high pressure, it should additionally consider the effect of pressure (e.g., Lei et al., 2016). Liquid-phase density ρ l takes into account the effects of the total pressure P, temperature T and CO 2 mass fraction in liquid phase x CO2 l . Both gas-phase density ρ g and CO 2 mass fraction in gas phase y CO2 g are functions of temperature T and both the partial pressures of CO 2 and H 2 O. The enthalpy of liquid-phase mixture can be written as where h CO2 is the enthalpy of CO 2 and is a function of the partial pressures of CO 2 and temperature. The heat of a solution of CO 2 in water expressed as h SOL is only related to temperature. The enthalpy of gas-phase mixture can be written as The void fraction α is a key parameter for two-phase flow, and can be expressed as For a heterogeneous mixture (i.e., u g ≠ u l ), the relation between the void fraction α and velocities is required. The commonly used model for heterogeneous calculations is the drift flux model developed primarily by Zuber and Findlay (1965). It is expressed as where the volumetric flux j is the volumetrically weighted velocity determined as j αu g + (1 − α)u l . The distribution parameter C 0 and drift velocity u d are affected by the local gas saturation, the equation for the void fraction can now be expressed as

Drift flux model
The parameters of C 0 and u d have significant effect on the void fraction α and further affect the pressure and temperature profiles. For different flow patterns (i.e., bubbly, slug, churn, annular), the values of C 0 ; u d vary as well (Hasan and Kabir, 2010). The parameters of C 0 and u d were determined by multiple different models in order to include the effects of flow pattern, pipe diameter and flow mixture in the drift flux model. Out of these models, only three most commonly used models were incorporated into the flow model: Rouhani and Axelsson. (1970), Hibiki and Ishii. (2003), and Shi et al. (2005).

Rouhani and Axelsson model
A comparison of void fraction correlations based on an unbiased data set showed that the Rouhani and Axelsson. (1970) model had a good fitting result (Woldesemayat and Ghajar, 2007). It is expressed as where σ is the surface tension between gas and liquid phase.

Hibiki and Ishii model
The pipe diameter has important effect on flow pattern. In a large diameter pipe, slug bubbles cannot be sustained. Due to the interfacial instability, slug bubbles in such flows will disintegrate to cap bubbles. For a large diameter pipe, Hibiki and Ishii. (2003) developed and recommended the parameters as Frontiers in Earth Science frontiersin.org

Shi et al. model
Based on experimental data in a 15-cm-diameter, 11-m-long plexiglass pipe, Shi et al. (2005) determined drift-flux parameters by using an optimization technique that minimizes the difference between experimental and model predictions. The model is written as where.
-C max is a maximum profile parameter between 1.0 and 1.5 and is set to 1.2 in this study.
-η is a parameter reflecting the effects of the flow status on the profile parameters, defined as β−B 1−B . B 2 C max − 1.0667 is the threshold parameter above which C 0 starts to drop below C max . β is calculated as β max (α, α|um| u sgf ). The gas superficial velocity is calculated as u sgf K u ( Frontiers in Earth Science frontiersin.org K 1.53 α ≤ a 1 1.53 where a 1 and a 2 are two transition points of gas fraction, suggested by Shi et al. (2005)

Flowchart of calculation
A rigorous numerical approach to solve Eqs 1, 2; Eqs 3, 4 includes generating a mesh of the wellbore, establishing the discrete equations, and solving them by the Newton-Raphson iterative method, which was implemented in the T2Well simulator (Pan and Oldenburg, 2014). However, this method has two drawbacks. The first is that the initial conditions for the two-phase flow are very difficult to determine and different initial conditions could result in wildly different results. The second is the lack of robustness of the method; iterative convergence failure is often encountered in the phase change region. Following the method proposed by Barelli et al. (1982), an improved method is demonstrated in Figure 5. This method is based on the fact that the temperature decreases and gas phase fraction increases from the downhole to the wellhead. There are two loops in the calculation, one for temperature, the other for gas phase fraction. For a given temperature decrease, the first loop aims  to find the change in the z value. The second loop aims to find the change in the gas-phase fraction.

Verification
To verify the model, a compared wellbore flow model is taken from Vasini et al. (2018). The downhole CO 2 mass fraction is about 3%. Figure 6 shows that the calculations are generally consistent with the measurements and simulation from T2Well-EWASG. The small temperature deviation is due to the reason that the model in this study has ignored the heat exchange with surrounding rock. Compared with T2Well-EWASG, the model is simple and easy to implement.

Model setup
The parameters and boundary conditions for wellbore models are listed in Table 2. Due to the lack of downhole fluid samples   Frontiers in Earth Science frontiersin.org 11 during the testing, the key parameter of the CO 2 mass fraction at the bottom of casing pipe is unknown. However, according to the previous mechanical descaling, the flash depth of both the three wells is above the casing bottom, which means that the fluid is in single-phase state at the casing bottom. Therefore, it can be calculated by the Henry's law that the CO 2 mass fractions at the bottom of the casing pipe for ZK304, ZK303 and ZK324 cannot be higher than 0.08%, 0.06% and 0.12%, respectively. To account for the uncertainty of CO 2 content, three different cases were simulated. In Case 1, the fluid is fully saturated with CO 2 , In Case 2, the amount of CO 2 is half the amount of the fully saturated case and Case 3 has no CO 2 present. Because of the large production rate and the short wellbore casing, the heat exchange with the surrounding rock is ignored in the model. The calculation is implemented from the downhole to the wellhead. The grid size for single-phase calculation is set to be 1.0 m, while for two-phase calculation it varies depending on the given temperature decrease in Figure 5. The value of 0.1 or less is suggested.

Results
The results from the aforementioned three drift flux models are compared. Figure 7 shows that there is a noteworthy difference  (Table 3). The CO 2 content has a certain effect on the pressure and temperature profiles; a low CO 2 content results in a lower temperature and pressure. The calculated wellhead pressures according to the Shi et al. model are respectively 3.2, 3.1 and 3.0 bar for the three specified scenarios of CO 2 content (Case 1, 2 and 3. The corresponding temperatures are 135, 134°C and 133°C, respectively.

FIGURE 13
Effect of production rate on performance.

Frontiers in Earth Science
frontiersin.org 14 ZK303 well is 0.47 m, 1.38 times bigger than the other two wellbores. Due to the low CO 2 content, the presence of CO 2 has small effect on the pressure and temperature profiles.
Similar to ZK324, the result from the Shi et al. model show best prediction capability for ZK304 as well (Figure 9). The corresponding average errors for pressure and temperature are 3.34% and 0.41%, respectively (Table 3). Due to the same reason as ZK303, the effect of CO 2 content on the profiles of pressure and temperature is ignorable.
Based on the optimal model (i.e., the best suitable DFM model for fitting), the distribution of another parameter in the wellbore is further analyzed. Figure 10 shows that the flash depth (i.e., the gas void fraction is larger than zero) for ZK324 is about 44.9 m without CO 2 , while it is 58.1 with 0.06% CO 2 (Case 2). The wellhead gas mass fraction is about 0.035 and the corresponding gas void fraction is about 0.75. For the case of 0.06% CO 2 content (Case 2), the CO 2 mass fraction in the gas phase reaches 0.38 at the flash depth and gradually decreases the next 30 m upwards due to the increase in steam content. The same trend with depth can be observed with the CO 2 mass fraction in the liquid phase. Compared with the measured result of CO 2 content of averaged 0.15%wt at the wellhead measured by Zhao et al. (1998), the calculation of 1.6% for Case 2 is higher. The reason for this difference is the fact that the measurement was conducted at the pressure of 1 bar, while the calculation was done with a much higher pressure of 3 bar. The decreased pressure results in water transformation from the liquid to gas phase and therefore a low CO 2 content was measured. Above the flash depth, the value of the gas void fraction rapidly increases by more than 0.3 and the churn flow pattern (i.e., gas void fraction usually large than 0.25) is dominant. Figure 11 demonstrates that the flash depth for ZK303 is about 112.5 m without CO 2 , while it is 119.5 with 0.03% CO 2 (Case 2). The wellhead gas mass fraction is about 0.075 and the corresponding gas void fraction is about 0.80. The CO 2 mass fraction in the gas phase is about 0.21 at the flash depth for the case of 0.03% CO 2 . The changes in the CO 2 mass fraction in both the liquid and gas phase are small above the depth of 100 m, which indicates that the more severe wellbore scaling mostly occurs between the flash depth and a depth of a 100 m. For the case of 0.03% CO 2 (Case 2) in the wellbore of ZK303, the calculated CO 2 content of 0.4% is close to the measurement of 0.29% by Zhao et al. (1998). Figure 12 indicates that the flash depth for ZK304 is about 97.7 m without CO 2 and 107.2 m with 0.04% CO 2 . The wellhead gas mass fraction is about 0.065 and the corresponding gas void fraction is about 0.84. The maximum CO 2 mass fraction in the gas phase for the case of 0.04% CO 2 reaches 0.27 at the flash depth. The CO 2 mass fractions in both the liquid and gas phase vary dramatically from the flash point down to a depth of 90 m. For the fluid containing half of the CO 2 content of the fully saturated case, the calculated wellhead CO 2 mass fraction in the gas phase is 0.6%, five times of the measurement of 0.12%. The reason has been explained above.
Comparing the results of the three wells, it can be seen that the proportion of liquid evaporation to steam for the relatively short wellbore length is small and therefore the wellhead temperature is high. The corresponding wellhead pressure is high due to the small pressure loss.
6 Discussions 6.1 Effect of production rate Based the optimal model, the production performance is analyzed by means of the production rate. It can be seen from Figure 13 that both the wellhead pressure and temperature increase and the gas phase mass fraction decreases with the production rate. With the consideration that the geothermal field will deplete with long-term development due to the decrease in reservoir pressure and wellbore scaling, the production rate is set to half of the testing rate. The results shows that the wellhead temperatures for ZK304, ZK303 and ZK324 decrease by respectively 6.1, 5.4°C and 3.1°C compared to the testing rate. The pressures decrease by 0.44, 0.36 and 0.27 bar, respectively. This is because the change of production rate affects the distribution of two phases and the pressure profile, ultimately affecting the flash process.

Effect of surrounding rock heat exchange
Giving the parameters k cem 1.4 W/(°C m), k T 2.1 W/(°C m), r 1 0.34 m, r 2 0.39 m, ρ r 2600 kg/m 3 and c r 1000 J/(°C kg), the calculated transient heat exchange coefficient is showed in Figure 14A. At the beginning of production, the value is about 26 W/(m 2°C ). For long-term production, the value decreases to about 1 W/(m 2°C ). Figure 14B demonstrated that the differences in the wellhead pressure, temperature and gas-phase mass fraction are small, within 0.008 bar, 0.12°C and 0.1%, respectively. This is because that the flow rate is large and flow path is short. Therefore, the heat loss can be ignored. However, we still suggest Frontiers in Earth Science frontiersin.org that discharge tests with different flow rates should be implemented to determine the effect of heat exchange.

Effect of calculation parameters
In order to analyze the effect of the calculation parameters such as the grid size and the temperature decrease, the additional two cases (Case 1: based case with dz=1.0, dt=0.1, Case 2: dz=0.5, dt=0.05, Case 3: dz=2, dt=0.2) are used for comparison ( Figure 15). The pressure and temperature differences are within 0.2 bar and 2.75°C, respectively. The maximal difference in gas mass fractions is about 1% at the wellhead. The difference in flash depths is small. This is because that it only determined by the value of dz. Therefore, a small grid size and temperature decrease should be used in the calculation, especially for the evaluation of gas mass fractions.
6.4 Application prospect of the two-phase wellbore flow model in geothermal field The two-phase wellbore flow model constructed by this paper has great application prospect. Firstly, it can be used to estimate the production performance parameters (such as wellhead temperature, pressure and dryness) which is very useful for new geothermal wells drilling. It can be coupled with reservoir simulators to evaluate the change in wellhead temperature, pressure and dryness as the geothermal development. Secondly, it is basis for calcite scaling evaluation which is very important for the injection of chemical antiscaling inhibitor and descaling. The wellbore calcite scaling includes four important coupled process: two-phase flow with phase transition, water-gas-mineral reaction, solute transport and wall adhesion ( Figure 16). Based on the these coupled models, the profile of the thickness of calcite scaling can be obtained, which is the future work. In the evaluation of calcite scaling, the two-phase wellbore flow would provide key parameters such as pressure, temperature, saturation and velocity for the following calculation of reactive transport. Thirdly, combination of the wellbore pressuretemperature measurement and two-phase model calculation, the CO 2 content in the reservoir can be inversely identified. Fourthly, it

FIGURE 15
Effect of calculation parameters based on ZK304.

FIGURE 16
Conceptual map for calcite scaling in wellbores.

Frontiers in Earth Science
frontiersin.org can be used to help design the separator at the wellhead and fluid (including gas and liquid phase) transportation pipeline.

Summary and conclusion
The Yangbajing geothermal field in Tibet is one of the typical high-temperature geothermal fields and is the first to be used for high-temperature power generation in China. During the development, numerous problems were encountered, such as resource depletion and wellbore scaling. In order to deeply understand these problems, this paper analyzes the two-phase flow in the typical production wells and the findings can be summarized as follows.
(1) Temperature and pressure measurements during the selfdischarge testing revealed that the temperature in the shallow geothermal reservoir is about 150°C and the pressure is around 7-8 bar at a depth of a 100 m. During the self-discharge testing, the wellhead pressure for all three wells was measured to be 2-3 bar. The temperature was 118°C-135°C at a depth of 10 m.
(2) A steady-state two-phase flow model for production wells is constructed based on the mass, energy and momentum equations. The difference between gas and liquid velocities (drift flux model) has a decisive impact on the temperature and pressure profiles. For a small diameter wellbore, the Shi et al. (2005) model has a better fitting result. For a large diameter wellbore, the Hibiki and Ishii. (2003) model shows a better performance. (3) According the optimal simulated results, the flash depths for ZK324, ZK303 and ZK304 without considering CO 2 are 44.9, 112.5 and 97.7 m, respectively. If the CO 2 content in fluid at the bottom of casing pipe is half of the saturation content, the flash depths would increase to 58.1, 119.5 and 107.2 m, respectively. The wellhead temperatures would be 134, 125°C and 128°C and the wellhead pressures would be 3.1, 2.3 and 2.6 bar, respectively. The wellhead gas-phase mass fractions for the three wells are between 3% and 8%. (4) Both wellhead temperature and pressure increase with the production rate, but the gas mass fraction decreases. When the production rate reduces to half the value of the testing rate, the wellhead temperature and pressure reduce by a range of 3°C-6°C and 0.25-0.45 bar, respectively.
The model in this study only considers a steady state and a transient two-phase flow model still needs to be developed in the future. Additionally, the wellbore model should be coupled with a reservoir model to evaluate the effects of reservoir pressure depletion on the production. Furthermore, the lack of CO 2 mass fraction measurements in the fluid samples at the downhole needs still needs addressing, as the CO 2 mass fraction is a key parameter for accurately evaluating the flash depth.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Funding
This work was jointly supported by the CNNC's centralized R & D project "Research on Key Technologies of geothermal exploration, development and utilization" and the Science and Technology Innovation Talents Project of Sichuan Province, China (No. 2022jdrc0027).

YX was employed by CNNP Kunhua Energy Development Co. Ltd
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.