PROCEEDINGS, Thirtieth Workshop on Geothermal Reservoir Engineering Stanford University, Stanford, California, January 31-February 2, 2005 SGP-TR-176

COUPLING THE HOLA WELLBORE SIMULATOR WITH TOUGH2 Ashish Bhat Daniel Swenson Shekhar Gosavi Kansas State University Mechanical and Nuclear Engineering Department Manhattan, KS-66506, USA. e-mail: [email protected]

ABSTRACT This paper describes a coupled model between the reservoir simulator TOUGH2 and wellbore simulator HOLA. The goal is to provide increased capability in TOUGH2 to better model flow in geothermal systems containing inclined wells with multiple feedzones. The model accounts for varying flowing bottomhole pressure and flow entering the wellbore from multiple feedzones at different depths, with the fluid in different thermodynamic states. This approach facilitates a more accurate simulation of behavior of the geothermal reservoir being exploited. The standard TOUGH2 input data is extended to support the new capability. Some sample problems are solved using coupled simulator and compared with the results obtained from the current deliverability model in TOUGH2. MOTIVATION An Enhanced Geothermal System (EGS) is being exploited at the East Flank area of Coso Geothermal field, California. It is an ideal testing site with high rock temperatures at depths less than 10,000 ft and a high degree of fracturing and tectonic stresses. However, some of the wells within this portion of the reservoir are relatively impermeable and exhibit a significant drawdown. Characterization of the flow near the wellbore, becomes even more important in such settings. Field data for the Coso wells also reveals the presence of two-phase flow and multiple feedzones. This leaves the standard coupled wellbore flow option in TOUGH2 unusable, since it is limited to wells with a single feedzone. Previous authors have developed coupled well and reservoir models. Murray and Gunn (1993) presented coupling between the TETRAD reservoir simulator

and WELLSIM wellbore simulator. They generate a series of wellbore tables using WELLSIM, which are then used for interpolation by TETRAD. Hadgu et al., (1995) describe a coupling of TOUGH2 and WFSA. The WFSA simulator is one of the three parts of the WELLSIM simulator package (Gunn and Freeston, 1991). However, discussion with the author indicated that development of this coupled package had stopped and that the software was not readily available. This paper describes coupling of TOUGH2 specifically with the HOLA simulator. We plan to extend this to the GWELL, GWNACL and HOLA family of wellbore simulation models, to form a complete package. All source code is available publicly for use by others. HOLA WELLBORE SIMULATOR HOLA is a multi-feedzone geothermal wellbore simulator for pure water (modified after Bjornsson, 1987). It can handle both single phase and two phase flows in a well, with variable grid spacing and wellbore radii. GWELL(Aunzo,1990) and GWNACL are modified versions of HOLA, that can handle H2O-CO2 and H2O-NaCl systems respectively. They were developed using FORTRAN language. HOLA reproduces the measured pressure and temperature profile in a flowing well and determines thermodynamic properties of water, relative flow rates at each feedzone for a given discharge condition at the wellhead. It has two approaches (Option 1 and Option 2) for wellbore flow simulation. Option 1 needs known discharge condition at the wellhead (pressure, temperature and enthalpy), in addition to flow rates and enthalpies of all but the last feedzone. The simulator proceeds from wellhead-to-bottomhole to calculate the flowing temperature and pressure profile along the well. In Option 2, the user specifies the required flowing wellhead pressure and

bottomhole pressure and for each feedzone, the productivity indices and thermodynamic properties of reservoir fluid. The simulator then proceeds bottomhole-to-wellhead to calculate the expected wellhead output (wellhead enthalpy, flowrate, temperature and phase composition) for the required wellhead pressure. Governing equations are represented by two sets of equations. Namely, Between the feedzone and At the feedzone. The important equations are given below. Flow between feedzones Between the feedzones the flow is represented by one-dimensional steady-state momentum, energy and mass flux balances. Mass Balance •

dm =0 dL

(1)

•

Where

m = mass flow rate, L = length of pipe

Momentum Balance The total pressure gradient is the sum of the friction gradient, acceleration gradient and potential gradient.

dP ⎡ dP ⎤ ⎡ dP ⎤ ⎡ dP ⎤ −⎢ ⎥ −⎢ ⎥ −⎢ ⎥ = 0 dL ⎣ dL ⎦ fri ⎣ dL ⎦ acc ⎣ dL ⎦ pot

(2)

Where,

⎡ dP ⎤ 2 ⎡ dP ⎤ ⎢⎣ dL ⎥⎦ = φ Flow ⎢⎣ dL ⎥⎦ fri LO

(3)

d (G u m ) ⎡ dP ⎤ ⎢ dL ⎥ = dL ⎣ ⎦ acc

(4)

⎡ dP ⎤ ⎢⎣ dL ⎥⎦ = p g sin θ pot

(5)

φ2

Flow

[

θ

]

• dEt d 2 =m hm + 0.5u m + g (L − D ) (7) dL dL

Where hm is enthalpy of mixture, D is measured depth till the current grid-node. An approximate solution for Qt when we have the term, ∩ t / rw given as,

2

>> 1 (Carslaw and Jaeger, 1959) is

⎡ ⎧4 ∩ t ⎫⎤ Qt ≈ 4τ π (Tw − Tr ) ⎢ln ⎨ 2 − 2η ⎬⎥ ⎢⎣ ⎩ rw ⎭⎥⎦

−1

(8)

Where, η is the Euler’s constant (= 0.577216...), τ is rock thermal conductivity, t is time, ∩ is rock thermal diffusivity, Tw and Tr are temperatures in the well and reservoir respectively. Above equation does not take into account, transient changes in temperature and additional heat losses due to convection, when the well is flowing. However the heat loss term in equation (6) is very small compared to total energy flux term. At the feedzone At the feedzone, the mass and energy of inflow (or outflow) are given, and then mass and energy balance are performed, to continue further along the well. Here assumptions made are, instantaneous mixing occurs and it occurs at the wellbore pressure.

r m& m

r r = m& w − m& f

(9)

Subscripts m, w, f represent mixture, well, feedzone respectively. Flow from the feedzone for Option 2 is calculated using Darcy’s law as follows,

is the two-phase multiplier.

G is the mass flux, um is average fluid velocity, well deviation angle from horizontal.

dEt is the sum of dL

discharges in the heat content of the fluid, kinetic and potential energy. It is expressed as,

Mass Balance

⎡ dP ⎤ ⎢⎣ dL ⎥⎦ is the pressure drop for a flowing singleLO phase liquid and

The total energy flux gradient,

is

Energy Balance

⎡ k ρ k ρ ⎤ ⎡ dP ⎤ q = k A ⎢ rl l + rv v ⎥ ⎢ ⎥ µ v ⎦ ⎣ dr ⎦ ⎣ µl

(10)

where k is intrinsic permeability, krl and krv is relative permeability for water and vapor. µ is the viscosity.

dEt ± Qt = 0 dL

(6) Energy Balance

r m& m H m

r r = m& w H w − m& f H f

(11)

Where H is fluid enthalpy. Positive flowrate at the wellhead or a feedzone indicates production, while negative flowrate indicates injection. The mass flow in the well can have two possible directions, upward (producing) and downward (injecting). Also the flow for feedzone can be towards the well(producing) or towards reservoir (injecting). HOLA takes all such six combinations into consideration. It can be noted here that simulator WFSA assumes the fluid that enters the well from a feedzone, flows upward only. Formulae for two-phase flow calculations are taken after Chisholm (1983). Gas and liquid phase velocities are needed in the evaluation of momentum flux and energy equations. They can only be obtained by empirical correlations. Two choices are provided to user, to calculate these velocities, namely Armand(1946) correlation and Orkiszewski(1983) correlation. Reader is referred to Aunzo et al. (1991) for further details on input-output format of HOLA and various correlations. TOUGH2 SIMULATOR

TOUGH2 is a numerical simulator for nonisothermal flow of multi-component, multiphase fluids in one, two or three dimensional porous and fractured media (Pruess et al., 1999). It is a member of MULCOM family of codes, written in standard FORTRAN77 language and widely used in geothermal applications. Hence just a brief description of it is being given here with emphasis on sink/source feature. Darcy’s law is used to describe the single phase and two-phase flow. For the EOS1 water module, primary thermodynamic variables are pressure and temperature for a single phase flow, while pressure, temperature and phase saturation for two-phase flow. Various relative permeability functions and capillary pressure function, as a function of phase saturation are provided. TOUGH2 solves fluid motion by space discretization using the “Integral Finite Difference” Method (Edwards, 1972; Narasimhan and Witherspoon, 1976). Time is also discretized in a fully implicit manner as a first order backward finite difference. As a result, a set of strongly coupled nonlinear algebraic equations are formed, which are rearranged in a residual form. These equations are then solved iteratively by Newton-Raphson method, with time-dependent primary thermodynamic variables of all grid blocks as unknowns to be determined by bringing residual within a specified error limit. A well is represented as a sink/source. Sink/Source is specified directly as a grid block that acts as a region where injection or production of fluid mass and/or heat from the reservoir occurs. Various options are available for specifying the injection/production. One

of them is also a deliverability model for production wells to evaluate well-output based on a fixed specified bottomhole pressure and productivity index. The mass production rate with this option for phase β from a grid block with phase pressure P β > Pwb is,

qβ

=

k rβ

µβ

ρ β ⋅PI ⋅ ( Pβ − Pwb )

(12)

Thus the total rate of production for mass component κ is,

qˆ κ =

∑β X βκ qβ

(13)

X βκ is the mass fraction of phase β of component κ. As opposed to HOLA, a negative flow rate in TOUGH2 indicates production and a positive flow means injection. COUPLING OF HOLA WITH TOUGH2

Typically geothermal wells operate at almost a constant wellhead pressure. The well-output (flow rate and enthalpy at wellhead) varies with time, which means well bore pressure gradient and/or bottomhole pressure change with time. Thus it is more accurate to simulate a flowing wellbore against a variable bottomhole pressure. This exact approach is followed in our work, summarized as follows, i) A separate input file for each well in the reservoir is read in TOUGH2. ii) Unlike the current coupled wellbore flow option in TOUGH2 with ‘F----‘ type of wells (in which interpolation into predetermined well-tables is performed), an explicit call is made to HOLA, at the start of each new timestep. It can be noted that HOLA has been converted into a subroutine in the coupled code. iii) Required reservoir parameters are supplied to HOLA as input. Choice of initial guess value for bottomhole pressure and pressure step for iteration is left to the user through the input file. iv) HOLA then iterates using ‘option 2’ to calculate flow rates and enthalpies at various feedzones and at the wellhead. v) A positive (negative) flow rate calculated at a feedzone in HOLA is supplied as a constant production(injection) rate to the corresponding source/sink element in TOUGH2, during regular timestep calculations of TOUGH2. vi) Enthalpies for the producing element are calculated inside TOUGH2. For injection, enthalpy is taken from by HOLA. This approach is comparable to the two basic options in TOUGH2 for specifying the source/sink, namely ‘MASS’ option and ‘COM’ option. ‘MASS’ option specifies constant mass production rate with enthalpies to be determined from conditions in producing block. While option ‘COM’

represents mass injection at a constant rate with specified enthalpy. vii) Procedure from step (ii) to (vi) is then repeated for the next time-step with updated values of reservoir parameters. Now the initial guess for the bottomhole pressure is the converged value at the previous timestep. Some of the features of the coupled code are, • There is no change in the TOUGH2 input file. A coupled simulation with HOLA, is indicated in the GENER block of TOUGH2 input file. User specifies up to a five character long string name that starts with either ‘H’ or ‘h’, for each well. e.g. if a TYPE and ITAB record were specified as ‘H349’ and ‘a’, the program expects an input file with the name, ‘H349a’ present in the folder. • Input file for each well, is in similar spirit of the standard HOLA input file for option 2. • The user can also specify the wellhead pressure as time dependent data. Triple linear interpolation is then performed on this data and an average of wellhead pressure at the starting and ending time of current time-step is used for the calculation. • User may not want to keep a well flowing through the complete time-span of a reservoir run. Hence shut-in and flowing option is provided to the user. • With a proper choice pressure-step for iterations and a reasonable error limit on wellhead pressure, HOLA computes very fast, Hence comparatively there is not much compromise on computational time. The coupling procedure was aimed to make minimal changes in TOUGH2 and HOLA. The TOUGH2 subroutines in which the changes were made are, INPUT, RFILE, CYCIT, QU and OUT. Mostly new code blocks were added rather than a change in the original code. The input file for HOLA is read in RFILE and care is taken to properly preserve the sink/source indices and total number of generation grid blocks. The call to HOLA is made in CYCIT before the iterations for TOUGH2 begin. One more case was added to sink/source subroutine, QU to handle the coupled simulation with HOLA. Appropriate terms for Newton-Raphson iterations are calculated here. A small block of code was added to OUT, to get a printed output of well flow-rates, enthalpies and bottomhole pressure. The “Main” program of the original HOLA code is abandoned and a new subroutine HOLA( ) is written. Its purpose is same as the original main program, but there is no user interactive feature involved. Also some other new routines and functions are written. Here is their brief description,

i) SUBROUTINE PrepHOLA( ): gets the required reservoir parameter values from TOUGH2 and supplies them to HOLA. ii) SUBROUTINE TabPwh( ): finds current wellhead pressure from time-dependent table data. iii) SUBROUTINE FindTime( ): locates the current time in a tabular time data. iv) SUBROUTINE PwhInter( ): calculates wellhead pressure at given time, by doing triple linear interpolation between the values for starting time and end time. The hard-coded simulation parameters in HOLA were all brought at one place, so that they are accessible to user and could be changed as per problem size. Choice of a pressure step for bottomhole pressure iteration is left to the user’s discretion. More study is required to improve the performance of HOLA (Hadgu and Bodvarsson, 1992). Some issues in HOLA that have been addressed are, i) In subroutine VINNA2, while calculating the mass flow rate for the last feedzone, an average of reservoir fluid parameters and wellbore fluid parameters such as density, saturation, viscosity was being taken. A scenario in which one fluid could be single phase and other two-phase was not handled, leading to erroneous results. As an alternative to that we followed a simpler approach of using only reservoir parameters when fluid is entering the well and for the case of injection we use only wellbore fluid parameters. ii) Some instances of un-initialized variables have been corrected. iii) For relative permeability calculations, HOLA follows a simple approach. They are assumed to vary linearly between 0 to 1 and are equal to the phase saturation. We extended the various relative permeability calculation options available in TOUGH2 to HOLA also. And then HOLA uses the same choice of calculation as chosen by user in TOUGH2 input file. iv) We encountered a case in which the reservoir pressure was same as iterative value of bottomhole pressure. It lead to a division by zero calculation putting the program into an infinite loop. This was taken care of. v) We are correcting a minor error regarding the way inclined wells are currently handled in HOLA. SAMPLE PROBLEM

Sample problem number 5 from the TOUGH2 user’s guide is chosen as an example. This problem was originally taken after Hadgu et al., 1995. In this example, a coupled simulation model is compared with the deliverability option in TOUGH2.

A well of inside diameter 0.2 m produces from a 500 m thick two-phase reservoir containing water at initial conditions of P = 60 bars, T = Tsat (P) = 275.5 ˚C, Sg = 0.1. Wellhead pressure is 7 bars and the feedzone depth is 1000 m. A 1-D radial grid was created by means of MESHMAKER input file provided in the TOUGH2 user’s guide. The wellblock radius is 100 m and the grid extends to 10,000 m. The well productivity index is 4.64e-11 m3. Simulation starts with a time-step of 1.e5 seconds and ends at a time of 1.e9 seconds (approx. 31.7 years) First a run with coupled simulation is performed. Then ‘H----’ type in GENER block is replaced by DELV type so as to run the problem on deliverability model, with a fixed bottomhole pressure of 57.12 bar. This was the bottomhole pressure at the end of first timestep in the coupled simulation. Results obtained are plotted as shown in Figure 1. In the figure, Q is flow rate and h is the flowing enthalpy. Both Q and h are plotted for the coupled simulation (HOLA) and the deliverability model (DELV). Pwb (HOLA) is the varying bottomhole pressure for coupled simulation.

1280 1270 1260

50

40

1250

Q(HOLA) Q(DELV) Pwb(HOLA) h(HOLA) h(DELV)

1240

Enthalpy (KJ/kg)

Flow rate(kg/s) or Pressure (bars)

1290

70

60

1230

30 1220

20 1.00E+05

1.00E+06

1.00E+07

1.00E+08

ACKNOWLEDGEMENT

The authors appreciate discussions with Teklu Hadgu at Sandia National Laboratories and Karsten Pruess at LBNL. This work is supported by the U.S. Department of Energy, Assistant Secretary for Energy Efficiency and Renewable Energy, under DOE Financial Assistance Award DE-FC0701ID14186. REFERENCES

Au, A., “TETRAD User Manual”, ADA International Consulting Inc., Calgary, Alberta, Canada, 1995. Aunzo, Z. P., “GWELL: A Multi-Component MultiFeedzone Geothermal Wellbore Simulator.”, M.S. Thesis, University of California at Berkley, Berkley, CA, USA, May, 1990. Aunzo, Z.P., Bjornsson, G., and Bodvarsson., G.S. , “Wellbore Models GWELL, GWNACL, and HOLA”, Lawrence Berkeley National Laboratory Report LBL-31428, Berkeley, CA, October 1991. Bjornsson, G., “A Multi-Feedzone, Geothermal Wellbore Simulator.”, Earth Sciences Division, Lawrence Berkeley Laboratory Report LBL-23546, Berkeley, USA, 1987.

1300

80

GWELL and GWNACL simulators to form a complete package of coupled code.

1210 1.00E+09

Time(sec)

Figure1 :Flow rate, enthalpy and bottomhole pressure plotted against a logarithmic time scale. These trends match the results given in the TOUGH2 users manual. The deliverability model shows a rapidly declining production rate while the coupled reservoir-wellbore system shows a long term production at a much higher rate than a deliverability model. This certainly emphasizes the necessity of a coupled-reservoir model. CONCLUSION

With the results obtained we believe that we should be able to predict more accurately the production at the Coso EGS site. Future plans involve a careful review of the HOLA source code and coupling of

Carslaw, H. S., and Jaeger, J. C., “Conduction of Heat in Solids”, 2nd ed., Oxford Univ. Press, London, 1959. Edwards, A.L., “TRUMP: A Computer Program for Transient and Steady State Temperature Distributions in Multidimensional Systems”, National Technical Information Service, National Bureau of Standards, Springfield, VA, 1972. Gunn, C., and Freeston, D., “Applicability of Geothermal Inflow Performance and Quadratic Drawdown Relationships to Wellbore Output Curve Prediction”, Geothermal Resources Council Transactions, Vol. 15, pp. 471-475, 1991. Hadgu, T., Bodvarsson, G.S., "Supplement to wellbore models GWELL, GWNACL, and HOLA User’s Guide," LBL report LBL-32907, 1992. Hadgu, T., Zimmerman, R.W. , and Bodvarsson., G. S., “Coupled Reservoir-Wellbore Simulation of Geothermal Reservoir Behavior,” Geothermics, 24(2), 145—166. LBL-36141, 1995. Murray, L., and Gunn, C., “Toward Integrating Geothermal Reservoir and Wellbore Simulation:

TETRAD and WELLSIM”, Proceedings of the 15th New Zealand Geothermal Workshop, Geothermal Institute, The University of Auckland, Auckland, New Zealand, 1993. Narasimhan, T.N., and Witherspoon, P.A., “An Integrated Finite Difference Method for Analyzing Fluid Flow in Porous Media”, Water Resources Research, Vol. 12(1), pp. 57-64, 1976. Orkiszewski, J., “Predicting Two-Phase Pressure Drops in Vertical Pipe”, J. Pet. Tech., pp. 829-838, 1967. Pruess, K., Oldenburg, C., and Moridis, G., “TOUGH2 User's Guide, Version 2.0”, Lawrence Berkeley National Laboratory Report LBNL-43134, Berkeley, CA, November 1999.