Geothermal Wells - Two Phase Flow Modeling

Geothermics 36 (2007) 243–264

Method for selecting casing diameters in wells producing low-enthalpy geothermal waters containing dissolved carbon dioxide Vassilios C. Kelessidis a,∗ , Grigorios I. Karydakis b , Nikolaos Andritsos c a

c

Mineral Resources Engineering Department, Technical University of Crete, Polytechnic City, 73100 Chania, Greece b Institute of Geological and Mineral Exploration, Mesogeion 70, 11527 Athens, Greece Department of Mechanical & Industrial Engineering, University of Thessaly, Pedion Areos, 383 34 Volos, Greece Received 18 April 2006; accepted 18 January 2007 Available online 12 March 2007

Abstract Most low-enthalpy geothermal waters contain dissolved gases (e.g., CO2 , H2 S, and CH4 ). In artesian geothermal wells, the absolute pressure of the water flowing towards the surface may drop below the bubble point of the dissolved gases, resulting in their gradual release and the appearance of two-phase flow. To optimize flow conditions we must keep frictional losses to a minimum and prevent undesirable flow regimes from occurring in the well. A mechanistic model has been developed for upward two-phase flow in vertical wells, based on existing correlations for the various flow regimes. Computations have been performed using data measured in wells at the Therma-Nigrita geothermal field, Greece. The methodology presented here allows us to study the effects of changes in well casing diameter on fluid production rate and flow stability within the well, parameters that have to be considered when designing geothermal wells for further exploitation and field development. © 2007 CNR. Published by Elsevier Ltd. All rights reserved. Keywords: Low-enthalpy geothermal wells; Fluid production; Two-phase flow; Well design

∗

Corresponding author. Tel.: +30 28210 37621; fax: +30 28210 37874. E-mail address: [email protected] (V.C. Kelessidis).

0375-6505/$30.00 © 2007 CNR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.geothermics.2007.01.003

244

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

Nomenclature A b,c dp/dz D Db Dc f F g G H ki KCO2 lE L Lw m ˙ p pb pCO2 pe pr p0 pi pLw pT ptotal Q ReT U UGS ULS xi X yi zi

pipe/casing cross-sectional area (m2 ) constants [Eqs. (20)–(23)] pressure gradient (atm/m) (1 atm = 0.101352 MPa) well string diameter (m) depth of the location of the bubble point (m) critical well string diameter [Eq. (4)] (m) friction factor mixture molar feed (inlet) rate (mol/s) acceleration of gravity (m/s2 ) gas molar rate (mol/s) height of gas or liquid column (m) equilibrium constant for species (i) ((1) H2 O; (2) CO2 ) Henry’s constant for CO2 (atm) distance to the location of the bubble point (m) liquid molar rate (mol/s) length of liquid column (m) mass flow rate (kg/s) absolute pressure (atm) bubble point pressure (atm) partial CO2 pressure (atm) pressure at the exit of section (i) (atm) reservoir pressure (atm) water vapor pressure (atm) total pressure loss for section (i) (atm) pressure loss for liquid only flow (atm) pressure loss for two-phase flow (atm) total pressure loss (atm) volumetric rate (m3 /h) Reynolds number for two-phase flow [Eq. (16)] average velocity (m/s) superficial gas velocity (m/s) superficial liquid velocity (m/s) molar fraction of species (i) in the feed (inlet) Lockhart–Martinelli parameter molar fraction of species (i) in the gas phase ((1) H2 O; (2) CO2 ) molar fraction of species (i) in the liquid phase ((1) H2 O; (2) CO2 )

Greek letters α gas void fraction μL liquid viscosity [kg/(m s)] νL liquid kinematic viscosity (m2 /s) ρ density (kg/m3 ) σ water-air interfacial tension (N/m) ΦLo parameter defined in Eq. (19)

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

245

Subscripts fr frictional fr-Lo frictional, liquid only fr-w frictional, water gr gravitational gr-w gravitational, water G gas L liquid T two-phase wh, wh2 wellhead

1. Introduction Low-temperature geothermal fluids (i.e., temperatures less than 90 ◦ C; Muffler and Cataldi, 1978) often contain significant amounts of dissolved gases at reservoir conditions. The main dissolved gas is usually CO2 , although other gases such as CH4 , H2 S and N2 may also be present. In Greece, for example, the majority of these low-enthalpy fluids contain CO2 at ratios of up to 7.8 g CO2 /kg H2 O (Andritsos et al., 1994). In all geothermal fields the presence of gases in reservoir fluids has to be considered when designing and implementing a drilling program. These gases also present significant challenges in the production and collection of the hot fluids and their transmission to the utilization plants. There are a number of important factors that have to be taken into account when designing, drilling and completing geothermal wells since the final objective is to achieve the maximum possible flow rate without any significant drop in the temperature of the produced fluids. According to Antics (1995) and Karydakis (2003) these factors are: (a) the selection of appropriate diameters and depths for surface and intermediate borehole casings; (b) good cementing of these casings to avoid inflow of lower temperature fluids into the wells; (c) the selection of appropriate diameters for the production casing, allowing maximum fluid flow rate at minimum frictional pressure loss and ensuring that undesirable two-phase flow patterns (slug, churn or annular flow) do not form in the production string; (d) the reduction of heat loss to the immediate surroundings; in low-temperature geothermal systems these well losses are generally insignificant. During production of low-enthalpy geothermal fluids, CO2 may be released and a two-phase flow may appear in the wellbore, in which the gas and liquid phases may assume different flow patterns. The typical flow patterns observed during vertical upward two-phase (gas–liquid) flow in pipes are shown in Fig. 1. For constant liquid flow rate and increasing gas rate, the flow patterns are bubble, slug, churn and annular flow (Taitel et al., 1980; Hewitt, 1982), although periodic flows (geysering) are also possible (Lu et al., 2005). Apart from the diameter and roughness of the casing, the frictional pressure losses also depend on the particular flow pattern affecting absolute pressures along the wellbore, which determine the amount of gas released from solution, further modifying the existing flow pattern. Hence, different patterns may develop along the well as the fluid ascends towards the surface (Szilas and Patsch, 1975; Garcia-Gutierrez et al., 2002).

246

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

Fig. 1. Schematic depiction of flow patterns for two-phase flow in pipes.

It is important to optimize the design of wells to be able to predict the flow patterns occurring in the boreholes in order to avoid fluid flow-related problems. One should remember that the behavior of the wells will definitely have an impact on the performance of the surface pipelines. The desirable flow patterns in a producing geothermal well are, in order of preference, single-phase (if possible), bubble, and dispersed bubble flow. Patterns that should be avoided, given in terms of increasing undesirability, are slug, churn and annular flow, mainly because fluid flow is more difficult to control. Not much has been published on modeling such flows in shallow geothermal wells, except for the studies by Tolivia (1972), Szilas and Patsch (1975), and Antics (1995). Gunn et al. (1992a,b) addressed the issues of calibrating and validating wellbore simulators for deep geothermal wells, while Garg et al. (2004) presented a new liquid hold-up correlation based on measurements for deeper wells in conjunction with a simulation code. Recently, Lu et al. (2006) discussed experimental and modeling results of transient two-phase flow in shallow geysering geothermal wells. Here we present (a) a model for the fluid mechanics of artesian low-enthalpy fluid production in vertical geothermal wells, (b) a comparison of our predictions with measurements in producing wells, and (c) a proposed methodology for optimizing the design of future drilling programs. The importance of such a methodology becomes evident if we consider how well construction costs dominate the economics of geothermal power generation (Combs et al., 1997); they can also represent a significant component of final electricity prices (Garg and Combs, 1997), and could account for 50% (Barbier, 2002) to 70% of the total cost of a geothermal project (Antics and Rosca, 2003).

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

247

Fig. 2. Schematic diagram of phase separation of geothermal fluids as they ascend to the surface.

2. Theory 2.1. Flow patterns in producing low-enthalpy artesian geothermal wells During production, the flow of low-enthalpy fluids in vertical well casing strings, typically 0.0762–0.2032 m (3–8 in.) in diameter, can be either single-phase or two-phase, depending on the prevailing conditions in the well and reservoir. Single-phase flow normally occurs in the lower part of the well, where, because of high pressures, the gases are in solution (Tolivia, 1972). As the geothermal fluid flows upwards towards the surface, the pressure in the fluid column decreases due to the smaller hydrostatic pressure and frictional losses. At some point in the well, the sum of the partial pressures of the dissolved gases may become equal to the absolute well pressure, at which point CO2 (and/or other gases) will start coming out of solution, gas bubbles will form, and bubble flow begins (Fig. 2). This flow pattern is characterized by discrete smalldiameter gas bubbles that move upwards in a zig-zag manner at a faster rate than the liquid. Further up the well, the absolute pressure decreases, resulting in the exsolution of more gas and its expansion, generating even larger bubbles. This increases bubble density in the mixture to the point where coalescence of the smaller bubbles results in the formation of larger bubbles (Taylor bubbles), causing the transition to slug flow (Taitel et al., 1980; Kelessidis and Dukler, 1989) (Fig. 1). Further up the well, more gas comes out of solution and gas expansion continues. The Taylor bubbles grow in length, increasing bubble velocity and total gas volumetric flow rate. As the Taylor bubbles ascend, the liquid falls between the pipe and the bubbles forming a film that penetrates deeply into the liquid slug following the Taylor bubbles, creating a gas–liquid mixture containing large amounts of gas; this results in the disintegration of the liquid slug and transition to churn flow (Fig. 1). Churn flow has been characterized as an entrance region phenomenon in vertical pipes (Taitel et al., 1980) and in vertical annuli (Kelessidis and Dukler, 1989), although there is still scientific debate about the existence of this particular flow regime (Jayanti and Hewitt, 1992; Chen and Brill, 1997). Reports from continuous monitoring and visual observations of two-phase (gas–liquid)

248

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

flow in vertical annuli indicate that, during churn flow conditions, the gas moves continuously upwards, lifting the liquid to a certain height; the liquid then falls, accumulates, bridges the pipe and is again lifted up by the gas (Kelessidis, 1986; Kelessidis and Dukler, 1989). This chaotic oscillatory motion of the liquid is the main characteristic of churn flow, which is expected to occur close to the location of the bubble point. Based on various reports and observations (e.g., Tolivia, 1972; Antics, 1995; Lu et al., 2005), it appears that churn flow is very likely to develop in low-enthalpy geothermal wells during production. The transition to annular flow occurs at very high gas flow rates (Fig. 1). Liquid ascends as a film covering the wall of the pipe, while gas flows upwards in the core carrying liquid droplets entrained from the liquid film. This flow pattern is not expected during the production of lowenthalpy fluids because the amount of gas in the gas–liquid mixture is never high enough for annular flow to exist. 2.2. Prediction of flow pattern transitions in pipes In most cases, flow pattern transitions are gradual as the liquid and gas phase flow rates change. When these transitions occur, the flow features of both patterns are often observed over a narrow range of flow rates (Kelessidis and Dukler, 1989). Such transitions are depicted in flow pattern maps that have as coordinates the superficial gas and liquid velocities, UGS and ULS , given by: UGS =

QG A

(1)

ULS =

QL A

(2)

and

where QG and QL are the gas and liquid phase volumetric rates and A is the pipe cross-sectional area. Examples of such maps are given in Section 3.3 below. Taitel et al. (1980) provided the most comprehensive models for flow pattern transitions for upward gas–liquid flow in pipes, while modifications for annulus geometry were presented by Kelessidis and Dukler (1989). The occurrence of a particular flow pattern depends on the void fraction, α, defined as: α=

UGS UG

(3)

where UG is the average cross-sectional gas velocity. For bubble flow to exist, the velocity of the bubbles must be smaller than the velocity of the Taylor bubbles. This gives a condition between the critical pipe diameter (Dc ), fluid and gas densities (ρL , ρG ), and liquid surface tension (σ) that is given by Taitel et al. (1980): 

ρL2 gDc2 σ(ρL − ρG )

1/4 = 4.36

(4)

where g is the acceleration of gravity. If the pipe diameter D is larger than Dc , then bubble flow will be observed; otherwise that pattern should not be expected.

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

249

At low liquid flow rates transition from bubble to slug flow takes place when (Taitel et al., 1980):   g(ρL − ρG )σ 1/4 1−α ULS = UGS − 1.53(1 − α)1.5 (5) α ρL2 This transition occurs when the void fraction, α, becomes equal to a specific value, with most researchers suggesting the value of 0.25 for flow for various conduits (Taitel et al., 1980; Kelessidis and Dukler, 1989). Thus, Eq. (5) becomes:   g(ρL − ρG )σ 1/4 ULS = 3.0UGS − 0.994 (6) ρL2 and this curve is denoted as (A) in a flow pattern map. At high liquid rates, turbulent forces break up the small gas bubbles, resulting in a finely dispersed bubble regime where the void fraction can exceed the value of 0.25 without observing a transition to slug flow. Taitel et al. (1980) proposed that this happens when:     D0.429 (σ/ρL )0.089 g(ρL − ρG ) 0.446 (7) ULS + UGS = 4.0 ρL νL0.072 where νL is the liquid kinematic viscosity. This equation is denoted as curve B in a flow pattern map, and cannot extend to values of the void fraction higher than the maximum packing of bubbles, which, for the case of cubic packing, occurs at a void fraction of 0.52. This leads to Eq. (8) below, derived from Eq. (5) for α = 0.52, and denoted as curve C in a flow pattern map:   g(ρL − ρG )σ 1/4 ULS = 0.9231UGS − 0.5088 (8) ρL2 It has been shown (Taitel et al., 1980) that churn flow will be observed at a distance lE from the location of the bubble point, if the gas and liquid superficial velocities satisfy Eq. (9):   ULS + UGS lE √ = 40.6 + 0.22 (9) D gD Eq. (9) is shown as curve D in a flow pattern map for a given value of lE /D. For the churn-to-annular flow transition, Taitel et al. (1980) proposed that annular flow cannot exist unless the velocity of the gas in the core is high enough to sustain the maximum size of the entrained liquid droplets, which is represented as curve E in a flow pattern map and is given by: 1/2

UGS ρG = 3.1 [g(ρL − ρG )σ]1/4

(10)

2.3. Estimation of pressure losses for two-phase flow in vertical pipes The overall pressure loss in non-horizontal pipes, after neglecting acceleration effects, is (Dukler and Taitel, 1986):     dp dp dp + (11) = dz dz fr dz gr where (dp/dz)fr and (dp/dz)gr are the frictional and hydrostatic pressure losses, respectively.

250

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

The pressure loss due to gravity is given by:   dp = g[ρL (1 − α) + ρG α] dz gr

(12)

The void fraction is computed from Eq. (3) and can also be related to the Lockhart–Martinelli parameter, X (Wallis, 1969): −0.378

α = (1 + X0.8 )

(13)

where X is defined as the square root of the ratio of the pressure loss in the pipe for liquid-only flow to the pressure loss in the pipe for gas-only flow (Lockhart and Martinelli, 1949). For two-phase bubble flow or dispersed bubble flow, the void fraction can be computed from Eq. (5), while the pressure loss is given by (Govier and Aziz, 1972):   dp 2fT ρL (ULS + UGS )2 (14) = dz fr D where fT is the two-phase friction factor, determined for turbulent flow from a Blasius-type equation (Wallis, 1969; Govier and Aziz, 1972): fT =

0.046 Re0.2 T

(15)

with the two-phase Reynolds number, ReT , defined as: ReT =

ρL D(ULS + UGS ) μL

(16)

where μL is the liquid viscosity. For slug flow, the frictional pressure loss is estimated as for bubble flow (Eq. (14)), with the void fraction computed by Eq. (3), but using as gas velocity, UG , the Taylor bubble velocity given by:  UG = 1.2(ULS + UGS ) + 0.35 gD (17) For churn flow, the frictional pressure loss is estimated as (Kern, 1975):     dp dp 2 = ΦLo dz fr dz fr-Lo

(18)

where ΦLo = cXb c=

14.2 (m ˙ L /1.64πD)0.1

b = 0.75

(19) (20) (21)

and m ˙ L is the liquid mass rate. The frictional pressure loss for liquid only (dp/dz)fr-Lo , is estimated using standard single-phase correlations (Govier and Aziz, 1972), such as Eqs. (14)–(16) but with liquid-only parameters (i.e., setting UGS = 0 m/s and ULS = UL ).

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

251

Wells producing low-enthalpy geothermal fluids normally do not present annular flow. For such a flow regime, the pressure loss is computed using the same equations as for churn flow but with the constants c, b given by (Kern, 1975):   D c = 4.8 − 0.3125 (22) 0.0254   D (23) b = 0.343 − 0.021 0.0254 2.4. Estimation of gas-phase concentration in a vertical well When the sum of the partial pressures of the non-condensable gases exceeds the fluid pressure at some point in the vertical pipe (i.e., vertical well), there will be a partial release of the dissolved gases and establishment of two-phase flow. Assuming that all of the dissolved gas is CO2 , the molar fraction of CO2 in the gas phase and the velocity of that phase can be determined following the procedure described below, which is based on vapor–liquid equilibrium considerations. Referring to Fig. 2, and for a geothermal fluid with a total molar flow rate F in the liquid state, containing two species, water (i = 1) and CO2 (i = 2), the exsolution of species 2 occurs somewhere between points A and B along the well, where the absolute pressures are pA and pB , respectively, and it holds that: pB < pCO2

(24)

where pCO2 is the partial pressure of CO2 . Total mass balance between feed (or inlet) point A and exit point B gives: F =G+L

(25)

with G, L the molar rate of gas and liquid, respectively, at point B. Mass balance for species (i) gives: xi F = yi G + zi L

(26)

with xi the molar fractions of species (i) in the liquid state (at point A), and yi , zi the molar fractions of species (i) in the gas and in the liquid phases at exit point B, respectively. The thermodynamic balance equation for species (i) is given by: ki =

yi zi

(27)

where ki is the equilibrium constant for species (i), which, for CO2 (i = 2), is given by: k2 =

KCO2 p

(28)

where KCO2 is Henry’s constant for CO2 and p is the absolute pressure of the fluid. The equilibrium constant for water (i = 1), with p0 the vapor pressure of water at the prevailing temperature, is: k1 =

p0 p

(29)

252

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

A combination of the above equations yields the gas molar rate: G=−

x1 (k1 − 1) + x2 (k2 − 1) F (k1 − 1)(k2 − 1)

(30)

while the molar fraction of water in the liquid phase at point B is given by: z1 =

x1 G(k1 − 1)/F + 1

(31)

and the molar fraction of water in the gas phase at point B by: y1 = z1 k1

Fig. 3. Flow diagram of the method used to estimate the location of the bubble point in a wellbore.

(32)

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

253

If the reservoir pressure, pr , and the wellhead pressure, pwh , are known, the depth at which the first gas exsolution occurs in the well (bubble or flashing point) can be determined by trial-and-error procedures. That depth can also be estimated from field measurements. If no measurements are available, we begin our computation of the bubble point by assuming that gas exsolution starts at a given depth, Db , where the pressure is equal to the bubble point pressure, pb . The pipe is subdivided into n sections (i.e., depth intervals) where conditions are assumed to be constant. From the top of the reservoir to Db , the total pressure loss for liquid-flow, pLw , is computed following the procedure described in Section 3.3. The gas superficial velocity, UGS , is calculated for each section (for which the prevailing flow pattern has been determined). The pressure loss for two phase flow, pT , is then computed using Eq. (11), calculating the frictional and gravitational contributions corresponding to the dominant flow pattern. Summation of all pressure loss contributions, pi , gives the total pressure loss, ptotal , up to this depth. From these computations, one can determine the pressure at the exit of the particular section, pe , which is equal to the pressure at the entrance of the next section (of the pipe) and hence the gas flow rate can be computed, which allows determination of the prevailing flow pattern in the section. The total pressure at the exit of the last section must equal the wellhead pressure. Where this holds true, the computation ends, otherwise the procedure is repeated (i.e., iterated). Schematically, the procedure is shown in Fig. 3. 3. Field data and computation of flow patterns along the well 3.1. The Therma-Nigrita geothermal field Data were collected from the low-enthalpy Therma-Nigrita geothermal field in northern Greece. The geology, and the data from four wells, are given in Fig. 4. The conglomerates and sandstones that host the geothermal reservoir rest on a strongly faulted metamorphic basement and are overlain by impermeable clay-sand sequences that act as caprock. The reservoir has been estimated to extend over an area of 12 km2 , with a thickness varying from 20 to 65 m; the top of the reservoir occurs between 70 m and 500 m depth. Measured reservoir fluid temperatures are in the 40–64 ◦ C range. The reservoir produces a two-phase CO2 –H2 O mixture under artesian conditions. Carbon dioxide concentrations in the produced geothermal fluid are in the 3–4 kg/tonne range. In most wells wellhead pressures (with the valves closed) are between 3 atm and 7 atm. 3.2. Pressure estimate for well TH-1 Data were collected in well TH-1 when the well was closed and during production; well characteristics are reported in Table 1. Reservoir temperatures were measured with an electricalresistance thermometer lowered into the well via an electric cable. Wellhead measurements (pressure and temperature) were made using the set-up shown in Fig. 5. The data presented in Table 2 were obtained under two different conditions: (a) the well was shut-in until the upper part was filled with gas, temperature equilibrium was attained, and the presence of water vapor had reached its minimum values; (b) the wellhead valve was opened, and CO2 was allowed to expand and discharge until the well was filled with liquid only. We will now describe the method used to compute the lengths of the liquid and gas columns. The reservoir top (point A in Fig. 6) is at 120 m depth. Because CO2 was trapped in the upper part of the well when the wellhead valve was closed, it is evident that at some depth (point B)

254

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

Fig. 4. Therma-Nigrita low-enthalpy geothermal field, northern Greece. Top: NNE-SSW geological cross-section. Middle: well characteristics. Bottom: map showing well locations and isotherms at reservoir level (in ◦ C); contour interval: 2 ◦ C.

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

255

Table 1 Characteristics of well TH-1 Parameter

Value

Total depth Thickness of fluid feed zone Length of 5-in. diameter casing (cemented) Length of casing (production string) Diameter of casing (production string) Perforated length Reservoir temperature (measured)a

135 m 15 m (120–135 m depth) 45 m (0–45 m depth) 135 m (0–135 m depth) 0.076 m (3 in.) 15 m (120–135 m depth) 59.4 ◦ C

a

Temperature was measured using a logging tool.

Fig. 5. Measurement set-up at the wellhead of well TH-1.

fluid pressure becomes equal to the partial pressure of CO2 (bubble point). That particular depth can be determined as follows. The partial pressure of CO2 , for the conditions of well TH-1, is determined from Henry’s law as: pCO2 = KCO2 x2 = (1640 atm)(0.00221 mol CO2 /mol H2 O) = 3.62 atm Table 2 Wellhead pressures and temperatures measured in TH-1 Before CO2 expansion Pressure Temperature After CO2 expansion Pressure (gage) Temperature

3.70 atm 8 ◦C 0.77 atm 24 ◦ C

(33)

256

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

Fig. 6. Completion and initial conditions in well TH-1. Bottom of the gas column is at 30 m depth.

where Henry’s constant, KCO2 , has been taken as 1640 atm (Ellis and Golding, 1963) at a temperature of 24 ◦ C (Table 2), the salt concentration in the water is 2.4 g/kg and the gas molar fraction x2 = 0.00221 (Table 3). This pressure is very close to the value measured at the TH-1 wellhead when the valve was closed (i.e., p = 3.70 atm).

Table 3 Wellhead conditions during production of well TH-1 Parameter

Value

Pressure Temperaturea Water density CO2 density Water volumetric flow rateb CO2 volumetric flow ratec Kinematic viscosity of water Water surface tension Volumetric concentration of non-condensable gasesd CO2 content Mass of dissolved CO2 at the exit Total dissolved solids (measured)e

1 atm 59.4 ◦ C 983.2 kg/m3 1.61 kg/m3 1.39 × 10−2 m3 /s = 50 m3 /h 4.06 × 10−2 m3 /s = 146.2 m3 /h 4.75 × 10−7 m2 /s 66.2 × 10−3 N/m 99.2% CO2 0.54 g CO2 /100 g H2 O = 0.00221 mol CO2 /mol H2 O 4.7 kg CO2 /m3 H2 O = 2.4 Nm3 /m3 2.4 g/L

a

Temperature was measured using a digital thermometer. Liquid volumetric flow rate was measured with a 4-in. turbine flow meter. c Gas volumetric flow rate was measured with a 4-in. orifice meter at the vapor outlet of the surface liquid–gas separator. d The gas content in the vapor phase was obtained after collecting the vapor phase in a gas bottle and analyzing for components in the laboratory the same day, and determining the amount of the gases, including CO2 , by gas chromatography. e Total dissolved solids were determined based on conductivity-meter measurements. b

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

257

When the wellhead valve is closed, the pressure at the top of the reservoir, pr , is essentially given by: pr = pA = ρL gHL + pwh

(34)

where pwh is the wellhead pressure, ρL the (average) density of the water and HL is the height of the liquid column. As mentioned earlier, when the wellhead valve is opened for a brief period, and before production starts, the entire column of the well consists of liquid only. The pressure at the top of the reservoir is then given by: pA = ρL g(HL + HG ) + pwh2

(35)

where pwh2 is the new measured wellhead pressure. The water density is computed for a temperature of 42 ◦ C, the average between the measured bottomhole (59.4 ◦ C) and wellhead (24 ◦ C) temperatures, as ρL = 992.2 kg/m3 . Combining Eqs. (34) and (35) with the measured values of pwh and pwh2 yields HL = 90 m, HG = 30 m and pA = 12.3 atm. These measurements allowed us to estimate the pressure at the top of the reservoir, a value to be utilized in calculations related to the production phase. 3.3. Two-phase flow during production from well TH-1 Mass flow-rate measurements of liquid and gas were made in well TH-1 during production using a surface separator, a turbine flowmeter for liquid and an orifice meter for gas, while recording the pressure and temperature. The well schematics for this condition are shown in Fig. 7 and the data collected are given in Table 3. The liquid volumetric flow rate, measured with the flowmeter, together with the fluid density estimated from the temperature, can give the liquid mass flow rate. Likewise, the gaseous volumetric flow rate, consisting mainly of CO2 (Table 3),

Fig. 7. Prevailing flow patterns within the borehole of well TH-1 during production. Well diameter: 0.076 m.

258

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

obtained with the orifice meter at the gas outlet of the surface liquid–gas separator, using the measured pressure and temperature, can give the CO2 mass flow rates. In the lower part of the well, where the pressure is greater than the partial pressure of CO2 (i.e., p > pCO2 ), there is only single-phase flow. The bubble point occurs at point B in the wellbore, where total pressure equals the bubble pressure, i.e., p = pb = pCO2 . Taking Henry’s constant as KCO2 = 3440 atm for the well temperature of 59.4 ◦ C and x2 = 0.00221 from Table 3, then pCO2 = 7.6 atm. The length over which single-phase (liquid) flow exists is obtained by calculating the pressure loss from the top of the reservoir to the bubble point (p)Lw : pLw = pr − pCO2 = 12.3 − 7.6 = 4.7 atm

(36)

This pressure loss is equated to the gravitational, pgr-w , and frictional pressure loss, pfr-w pLw = pgr-w + pfr-w = ρL gLw +

2fρL UL2 Lw D

(37)

where Lw is the required length of liquid column. Using ρL = 983.2 kg/m3 (for 59.4 ◦ C), and taking the value of the friction factor (from the Moody diagram given in Govier and Aziz, 1972) as f = 0.006, for the case of turbulent flow of water in a steel pipe with roughness 0.15 mm and D = 0.076 m, the length for single-phase flow is estimated as Lw = 48.75 m. At depths shallower than (120–48.75) = 71.25 m, bubble flow should occur, assuming that the conditions for a bubble flow regime are satisfied. Based on Eq. (4) and using data from Table 3, Dc = 0.054 m; since D > Dc , bubble flow will exist above 71.25 m depth. Using the equations presented above, the gas exsolution rate, the mole fractions in the gas and the liquid phase, the superficial gas velocity, UGS , can be calculated for every point along the well. The liquid superficial velocity is computed from the measured rate at the wellhead, and as the temperature in the well does not significantly change during production, it remains constant at ULS = 3.11 m/s. Part of these computations is shown in Table 4 and the flow patterns prevailing in the well are indicated in Fig. 7. Bubble flow exists from 71.25 m (point B in Fig. 7) to 36.45 m depth (point C), where pressure is 4.31 atm. For this depth lE = (71.25–36.45) m and lE /D = 457. According to Eq. (9), for this ratio, ULS + UGS = 9.54 m/s, which is much higher than the value Table 4 Superficial gas velocity, gas content and flow pattern at different depths in well TH-1 Depth (m)

Pressure (atm)

Gas molar rate (mol/s)

zCO2 (×10−3 )

yCO2

UGS (m/s)

Flow pattern (point in Fig. 7)

71.25 65.38 55.14 44.34 36.45 31.61 20.49 10.64 5.54 0.00

7.60 7.00 6.00 5.00 4.31 3.88 3.00 2.25 1.80 1.24

– 0.145 0.378 0.613 0.777 0.880 1.097 1.293 1.422 1.611

2.23 2.02 1.73 1.43 1.23 1.11 0.85 0.63 0.49 0.33

– 0.983 0.980 0.976 0.972 0.969 0.959 0.946 0.932 0.902

– 0.125 0.380 0.738 1.086 1.367 2.203 3.462 4.756 7.824

Bubble flow (B) Bubble flow Bubble flow Bubble flow Slug flow (C) Dispersed bubble flow (D) Dispersed bubble flow Slug flow (E) Slug flow Slug flow (F)

Well diameter = 0.076 m. ULS = 3.114 m/s. UGS : superficial gas velocity; yCO2 : molar CO2 fraction in the gas phase; zCO2 : molar CO2 fraction in the liquid phase.

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

259

Fig. 8. Flow pattern map for production in well TH-1 assuming a well diameter of 0.076 m. Wellhead pressure: 1.24 atm; lE /D = 450 (see text for additional details). Curve A represents the points for transition from bubble to slug flow; curve B, from bubble to dispersed bubble flow; curve D from slug to churn flow; curve E, from churn to annular flow. Curve C corresponds to the points for which bubble flow can exist for the maximum void fraction value of 0.52. Curve F shows the actual flow patterns occurring in the TH-1 wellbore, starting from the bubble point (point G) to the wellhead (point H).

given in Table 4 for 36.45 m depth [i.e., 4.2 m/s = (3.114 + 1.086) m/s]. Hence, the flow pattern cannot be churn flow; it is in fact slug flow. At 31.61 m depth (point D), where pressure is 3.88 atm, the pattern changes to dispersed bubble flow, which persists to a depth of 10.64 m (point E), where the pressure is 2.25 atm and the regime changes again to slug flow. At the top of the well (i.e., wellhead; point F), we calculate a pressure of 1.24 atm, a superficial gas velocity of 7.8 m/s and slug flow. In reality, the wellhead pressure should have been 1.0 atm. The difference is attributed to the assumptions made in deriving the full model; however, the discrepancy is small and the results can be considered to be within engineering accuracy. The corresponding flow pattern map showing transition curves and the actual flow regimes in the TH-1 borehole during production is given in Fig. 8. Curve E is the transition curve most affected by pressure at a particular point, while curves A–D are not greatly affected by pressure; hence, the results shown in Fig. 8 can be considered a good representation of a flow pattern map at points lower in the well. Curve F represents actual computed values of the pair of superficial gas and liquid velocities, UGS and ULS , respectively, from the position of the bubble point, point G, to the wellhead, point H.1 Hence, above the bubble point, the flow patterns occurring in this well are bubble, slug, dispersed bubble and slug flow, as described before (Fig. 7). In well TH-1, at the wellhead, the distance (i.e., ratio) lE /D needed to develop churn flow is 71.25/0.076 = 937.5. Curve D shown in Fig. 8 is for lE /D = 450. Hence, churn flow should not occur in TH-1. 1 In Fig. 8, point G seems to correspond to U = 3.11 m/s and U LS GS = 0.01 m/s mainly because smaller log cycles are not given for clarity. Actually, point G corresponds to the tiniest gas superficial velocity, which in the representation of Fig. 8, starts at 0.01 m/s.

260

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

No direct flow measurements were made in well TH-1, but a velocity log was run in nearby well TH-8, which has similar characteristics. The methodology just described allowed us to predict bubble point at 83 m depth. The velocity log showed an increase in velocity around 85–90 m, indicating the location of the bubble point, which is in reasonable agreement with our predictions. 4. Estimation of well diameter for future drilling activity The theoretical model and methodology presented above allows us to predict the behavior of low-enthalpy wells during production. This information should prove useful when designing new wells and estimating their production capacities since it can be used to determine the sensitivity of the wells to certain design and production parameters. For example, the analysis of field data within the framework of the developed methodology may indicate unwelcome changes in the flow pattern along the wellbore and at the wellhead, as well as the likely occurrence of a fairly unstable flow pattern (i.e., slug flow). The diameter of future wells would then be duly increased so as to incur smaller pressure losses, and avoid such unfavorable fluid flow conditions. In terms of flow pattern stability, and based on the behavior of the two-phase mixtures, the most desirable flow pattern is dispersed bubble flow because the two-phase fluid forms a homogeneous mixture and the undesirable effects of the discrete phases, such as periodic or chaotic variation of gas and liquid flows and excessive pressure losses, may occur rarely. The methodology described here can be used to study the effects of well diameter, one of the most important parameters in the design of a geothermal borehole. An analysis was performed using data from well TH-1, assuming that the liquid and gas flow rates, and the gas concentration, remain the same in all cases considered. The results for a well diameter of 0.06 m (2.36 in.) are reported in Table 5 and the flow patterns occurring along the well from the bubble point (point G) to the wellhead (point H) are represented by curve F in Fig. 9. They show that, by decreasing well diameter and maintaining the same production rates, the calculated pressure is 1.01 atm at 22.30 m depth. Under the assumed conditions, dispersed bubble flow occurs along most of the length of the wellbore. The bubble point is at 72.54 m depth; there is dispersed flow up to a depth of 27.19 m, slug flow up to 22.39 m depth, and annular flow up to 22.30 m. No churn flow should develop since lE /D at the wellhead is 937.5 (see Section 3.3). Table 5 Superficial gas velocity, gas content and flow pattern at different depths in well TH-1 Depth (m)

Pressure (atm)

Gas molar rate (mol/s)

zCO2 (×10−3 )

yCO2

UGS (m/s)

Flow pattern

72.54 67.50 58.95 50.22 41.42 32.85 27.19 25.47 22.39 22.30

7.60 7.00 6.00 5.00 4.00 3.00 2.27 2.00 1.37 1.01

– 0.145 0.378 0.613 0.851 1.097 1.288 1.363 1.562 1.716

2.23 2.020 1.730 1.430 1.140 0.850 0.630 0.550 0.370 0.260

– 0.983 0.980 0.976 0.970 0.959 0.946 0.939 0.911 0.878

– 0.204 0.609 1.185 2.057 3.535 5.482 6.585 11.017 16.580

Dispersed bubble flow Dispersed bubble flow Dispersed bubble flow Dispersed bubble flow Dispersed bubble flow Dispersed bubble flow Slug flow Slug flow Annular flow Annular flow

Well diameter = 0.06 m. ULS = 4.993 m/s. UGS : superficial gas velocity; yCO2 : molar CO2 fraction in the gas phase; zCO2 : molar CO2 fraction in the liquid phase.

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

261

Fig. 9. Flow pattern map for production in well TH-1 assuming a well diameter of 0.060 m. Wellhead pressure: 1.01 atm; lE /D = 450. See text and Fig. 8 for further details.

The results also indicate that assumption of the same mass rates for gas and liquid will not hold under the conditions just described because the decrease in pressure is too large. The maximum possible pressure loss of 11.3 atm, i.e. (12.3–1.0 atm), occurs over a depth of 50.24 m. The geothermal system will still produce geothermal fluids because the reservoir pressure is sufficient to keep the wells flowing. The wells will self-adjust by lowering the fluid production rates so to decrease the superficial gas and liquid velocities, thus reducing pressure losses in the wellbore. We can estimate the production rates under these self-adjusting conditions by assuming different values in the calculations so that the computed wellhead pressure is equal to 1 atm. Such a computation yielded a volumetric water production rate of 28.1 m3 /h, a decrease of about 44%. Computations assuming larger well diameters, but keeping the same water and gas flow rates, were also performed. The results are shown in Table 6 and Fig. 10 for a diameter of 0.127 m (5 in.). For such a well, only bubble and slug flow are predicted to occur in the borehole. The calculated Table 6 Superficial gas velocity, gas content and flow pattern at different depths in well TH-1 Depth (m)

Pressure (atm)

Gas molar rate (mol/s)

zCO2 (×10−3 )

yCO2

UGS (m/s)

Flow pattern

70.68 64.32 53.20 41.29 30.06 14.61 0.00

7.60 7.00 6.00 5.00 4.14 3.00 2.06

– 0.145 0.378 0.613 0.818 1.097 1.345

2.230 2.200 1.730 1.430 1.182 0.850 0.570

– 0.983 0.980 0.976 0.971 0.959 0.941

– 0.045 0.156 0.264 0.426 0.789 1.409

Bubble flow Bubble flow Bubble flow Bubble flow Slug flow Slug flow Slug flow

Well diameter = 0.127 m. ULS = 1.114 m/s. UGS : superficial gas velocity; yCO2 : molar CO2 fraction in the gas phase; zCO2 : molar CO2 fraction in the liquid phase.

262

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

Fig. 10. Flow pattern map for production in well TH-1 assuming a well diameter of 0.127 m. Wellhead pressure: 2.06 atm; lE /D = 450. See text and Fig. 8 for further details.

wellhead pressure was 2.06 atm, indicating that the system can produce more fluid than in the case of well TH-1, whose diameter is smaller (0.076 m). The calculations for the 0.127 m diameter well, assuming the gas-to-water ratio to be constant at the measured value of 0.54 g CO2 /100 g H2 O, a water density of 983.2 kg/m3 , and a wellhead pressure of 1.0 atm, give a water production rate of 209 m3 /h. In other words, by increasing the diameter of the well from 0.076 m to 0.127 m, water production could be increased from 50 m3 /h (Table 3) to 209 m3 /h (i.e, an increase of 318%). The effect of well diameter on water and CO2 production rates can be estimated using the methodology suggested above, keeping the gas-to-water mass ratio constant and equal to the

Fig. 11. Calculated water flow rate as a function of well diameter assuming the reservoir conditions of well TH-1 (e.g., gas-to-water ratio of 0.54 g CO2 /100 g H2 O, and a reservoir pressure of 12.3 atm).

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

263

value measured in well TH-1 (Table 3) and assuming a wellhead pressure of 1.0 atm. The results of one such computation are shown in Fig. 11, where the volumetric rate of water is shown as a function of the well diameter, keeping all other conditions and parameters constant and equal to those of well TH-1. One can observe that the wells, and thus the geothermal field, would be able to produce significant amounts of geothermal fluid if the casing diameters were increased. This holds provided that (1) the temperature remains the same, (2) isothermal conditions in the wells (and reservoir) prevail as assumed when deriving the data of the system, (3) the gas-to-water ratio does not change with production rate, which has been observed in practice, and (4) the production capacity of the geothermal reservoir permits it (this would have to be determined by carrying out well tests and modeling studies).

5. Conclusions A model of vertical geothermal wells has been developed that allows us to determine prevailing two-phase flow parameters during the production of low-enthalpy fluids that contain dissolved carbon dioxide. A systematic analysis of such wells has been performed. Similar methodologies for CO2 -containing, low-enthalpy geothermal fluids cannot be found in the published literature. Our model uses two-phase gas–liquid relationships to predict flow pattern transitions in the borehole. The point where bubble flow is first observed can be estimated from thermodynamic equilibrium data and from measured gas-to-liquid ratios. Otherwise, the model can be used to iteratively compute the location of the bubble point in the well. In this case, the pressure loss is calculated using single- and two-phase flow relationships corresponding to the prevailing flow patterns, while the gas concentration is computed from thermodynamic equilibrium data. Field measurements are reported from a well in northern Greece that produces low-enthalpy, CO2 -rich geothermal fluids. These data are utilized to estimate the pressure at the top of the reservoir. The pressures in the wellbore, the flow patterns occurring at different depths, and the associated pressure losses were also computed. The model shows that the conditions present in the studied well are suitable for bubble, dispersed bubble and slug flow, but not for the development of churn or annular flow. The suggested methodology allows us to study the effects of well diameter changes on the fluid production characteristics of low-enthalpy geothermal wells. This particular approach provides significant data that are not commonly utilized, but which can be used in the design of future wells and in the development strategy for a given geothermal area. As expected, larger diameter wells tend to produce greater volumes of fluids. However, when designing or sizing a well one should be aware of the flow patterns that could develop along the wellbore so that an optimum diameter can be chosen to achieve better flow stability during the production of two-phase fluids.

Acknowledgments The first author, V.C. Kelessidis, would like to dedicate this work to his Ph.D. advisor, the late Prof. A.E. Dukler, for guiding him to the wonderful world of two-phase flow. Part of the numerical code was developed by Mr. Y. Aspirtakis. The authors would also like to thank the anonymous reviewers and the editors of the journal for their valuable suggestions.

264

V.C. Kelessidis et al. / Geothermics 36 (2007) 243–264

References Andritsos, N., Karabelas, A.J., Emannouel, Ch., Karydakis, G., 1994. Characterization of fluids from low-enthalpy geothermal fields in Greece. In: Proceedings of the Communications of International Symposium Geothermics 94 in Europe, February. BRGM, Orl´eans, France, pp. 173–180. Antics, M., 1995. Modeling two phase flow in low temperature geothermal wells. In: Proceedings of the World Geothermal Congress. Florence, Italy, pp. 1905–1910. Antics, M., Rosca, M., 2003. Geothermal development in Romania. Geothermics 32, 361–370. Barbier, E., 2002. Geothermal energy technology and status: an overview. Renew. Sust. Energ. Rev. 6, 3–65. Chen, X.T., Brill, J.P., 1997. Slug to churn transition in upward vertical two-phase flow. Chem. Eng. Sci. 52, 4269–4272. Combs, J., Garg, S.K., Prichett, J.W., 1997. Geothermal slim holes for small off-grid power projects. Renew. Energ. 10, 389–402. Dukler, A.E., Taitel, Y., 1986. In: Hewitt, G.F., Delhaye, J.M., Zuber, N. (Eds.), Flow Pattern Transitions in Gas-Liquid Systems. Measurements and Modeling. Advances in Multiphase Flow, vol. 2. McGraw-Hill, New York, NY, USA, pp. 1–94. Ellis, A.J., Golding, R.M., 1963. The solubility of CO2 above 100 ◦ C in water and in sodium chloride solutions. Am. J. Sci. 261, 47–60. Garcia-Gutierrez, A., Espinosa-Paredes, G., Hernandez-Ramirez, I., 2002. Study on the flow production characteristics of deep geothermal wells. Geothermics 31, 141–167. Garg, S.K., Combs, J., 1997. Use of slim holes with liquid feedzones for geothermal reservoir assessment. Geothermics 26, 153–178. Garg, S.K., Pritchett, W., Alexander, J.H., 2004. A new liquid hold-up correlation for geothermal wells. Geothermics 33, 795–817. Govier, G.W., Aziz, K., 1972. The Flow of Complex Mixtures in Pipes. Krieger Pub. Co., New York, NY, USA, p. 792. Gunn, C.I.M., Freeston, D.H., Hadgu, T., 1992a. Principles for wellbore simulator validation and calibration using matching analysis-I. Analytical techniques. Geothermics 21, 341–361. Gunn, C.I.M., Freeston, D.H., Hadgu, T., 1992b. Principles for wellbore simulator validation and calibration using matching analysis-II. Case Study—Well Rotokawa 5, New Zealand. Geothermics 21, 363–376. Hewitt, G.F., 1982. Flow regimes. In: Hetsroni, G. (Ed.), Handbook of Multiphase Systems. Hemisphere, Washington, DC, USA, p. 1024. Jayanti, S., Hewitt, G.F., 1992. Prediction of the slug-to-churn flow transition in vertical two-phase flow. J. Multiphas. Flow 18, 847–860. Karydakis, G.I., 2003. Geothermal fields of low enthalpy in Northern Greece: drilling engineering, reservoir engineering, and two-phase flow of geothermal fluids. Ph.D. Dissertation. Democritus University of Thrace, Xanthi, Greece (In Greek), 347 pp. Kelessidis, V.C., 1986. Vertical upward gas-liquid flow in concentric and eccentric annuli. Ph.D. Dissertation. Univ. of Houston, Texas, USA, 314 pp. Kelessidis, V.C., Dukler, A.E., 1989. Modeling flow pattern transitions for upward gas-liquid flow in vertical concentric and eccentric annuli. Int. J. Multiphas. Flow 15, 173–191. Kern, R., 1975. Piping design for two-phase flow. Chem. Eng. 23, 145–151. Lockhart, L.W., Martinelli, R.C., 1949. Proposed correlation of data for isothermal two-phase two component flow in pipes. Chem. Eng. Prog. 45, 39–48. Lu, X., Watson, A., Gorin, A.V., Deans, J., 2005. Measurements in a low temperature CO2 -driven geysering well, viewed in relation to natural geysers. Geothermics 34, 389–410. Lu, X., Watson, A., Gorin, A.V., Deans, J., 2006. Experimental investigation and numerical modeling of transient two-phase flow in a geysering well. Geothermics 35, 409–427. Muffler, P., Cataldi, R., 1978. Methods for regional assessment of geothermal resources. Geothermics 7, 53–89. Szilas, A.P., Patsch, F., 1975. Flow in geothermal hot water wells. Geothermics 4, 79–88. Taitel, Y., Barnea, D., Dukler, A.E., 1980. Modeling flow pattern transitions for steady upward gas-liquid flow in vertical tubes. AIChE J. 26, 345–354. Tolivia, E., 1972. Flow in geothermal wells (an analytical study). Geothermics 1, 141–145. Wallis, G.B., 1969. One-Dimensional Two-Phase Flow. McGraw-Hill, New York, p. 409.