- Email: [email protected]
Modeling of coupled wellbore-reservoir flow in steam-like supercritical geothermal systems
Geothermics 86 (2020) 101793
Contents lists available at ScienceDirect
Geothermics journal homepage: www.elsevier.com/locate/geothermics
Modeling of coupled wellbore-reservoir flow in steam-like supercritical geothermal systems
T
Alfredo Battistellia,*, Stefan Finsterleb, Marica Marcolinic, Lehua Pand a
Independent Consultant, 60033, Chiaravalle, AN, Italy Finsterle GeoConsulting, Kensington, CA, 94708, USA c PRCFA Dept., Saipem SpA, Via Toniolo 1, 61032 Fano, PU, Italy d Energy Geosciences Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA, 94720, USA b
A R T I C LE I N FO
A B S T R A C T
Keywords: Coupled wellbore-reservoir flow Geothermal reservoirs iTOUGH2 Supercritical conditions TOUGH2 T2Well
The capabilities to efficiently model the transport of mass and heat in geothermal reservoirs at supercritical conditions are necessary to support the efforts to explore and utilize supercritical hydrothermal systems, as well as to extend the utilization of existing geothermal fields at greater depths, where supercritical temperatures might be encountered. Similar capabilities are needed to model the flow in wells producing supercritical fluids, possibly in a coupled way with the reservoir flow. As a step in this direction, the new EOS2H module for subcritical and steam-like super-critical H2O-CO2 mixtures has been linked to T2Well, the extension of the TOUGH2 numerical reservoir simulator developed for the transient modeling of fully coupled wellbore-reservoir flow. The resulting code is verified against simulations performed with other supercritical reservoir simulators and steadystate subcritical and supercritical wellbore flow simulations. It is then used to perform coupled wellbore-reservoir flow simulations at the well-sector scale relevant to the conditions found in well IDDP-1, Krafla, Iceland.
1. Introduction Supercritical fluids exist near magmatic heat sources, and supercritical conditions are often encountered beneath presently exploited high-temperature geothermal reservoirs. These supercritical ‘roots’ are often neglected when building 3D numerical models of exploited geothermal systems. One example is given by the 3D numerical model of the Menengai geothermal systems (Rift Valley, Kenya) built using TOUGH2-EOS2 (Pruess et al., 1999), which was artificially cut in the up-flow area at about 2000 m below ground level to avoid supercritical temperatures beneath the presently exploited two-phase steam dominated reservoir (Montegrossi et al., 2015). On the other hand, O’Sullivan et al. (2015) presented natural-state simulations at Menengai using a deeper model made possible by the supercritical version of the AUTOUGH simulator (Croucher and O’Sullivan, 2008). Supercritical Equations of State (EOS) for the TOUGH2 and iTOUGH2 (Finsterle et al., 2017) numerical reservoir simulators have been developed in the past considering pure water by Cox and Pruess (1990); Kissling (1995); Brikowski (2001), and more recently by Croucher and O’Sullivan (2008) and Magnusdottir and Finsterle (2015). Developments for mixtures have been presented by Kissling (2005) for
binary mixtures of H2O-NaCl up to 620 °C, and by McKibbin and McNabb (1999) for ternary mixtures of H2O-NaCl-CO2. O’Sullivan et al. (2016) developed an AUTOUGH-based simulator for water and air mixtures in order to incorporate shallow unsaturated zones. Modeling of coupled wellbore-reservoir flow in supercritical conditions is also required to interpret well testing operations, and in general to model wellbore flow effects on the performances of wells discharging from supercritical steam-like reservoirs. Morin et al. (2016) presented steady-state and transient simulations of wellbore flow inspired by the conditions expected for well Venelle 2 in the Larderello field, Italy, which was deepened within the DESCRAMBLED EU sponsored project aimed at experimenting with novel drilling technology and at investigating supercritical conditions in continental Europe. They used the commercial multiphase pipe flow simulator LedaFlow, employed also to simulate the flow in oil and gas production wells. Simulations performed included the heat exchange with surrounding rocks, assuming fixed bottomhole pressure and temperature (P&T) conditions. By simulating start-up and well shut-in operations, they observed that pressures, temperatures and flow rates are considerably affected by the thermal state of the well. As a first step in modeling coupled wellbore-reservoir flow in
⁎
Corresponding author at: AMBEN Dept., Saipem SpA, Via Toniolo 1, 61032 Fano, PU, Italy. E-mail addresses: [email protected] (A. Battistelli), stefan@finsterle-geoconsulting.com (S. Finsterle), [email protected] (M. Marcolini), [email protected] (L. Pan). https://doi.org/10.1016/j.geothermics.2019.101793 Received 25 June 2019; Received in revised form 22 November 2019; Accepted 18 December 2019 0375-6505/ © 2019 Published by Elsevier Ltd.
Geothermics 86 (2020) 101793
A. Battistelli, et al.
can also be used up to 140 MPa and 365 °C with minor errors (Battistelli et al., 2017). Density of pure CO2 is computed using the Peng and Robinson (1976) cubic EOS (PR EOS), while its dynamic viscosity is estimated using the so-called LBC method of Lohrentz et al. (1964). Pure CO2 enthalpy is computed by evaluating the ideal gas enthalpy with the Shomate equation (Chase, 1998), which is valid between 298 and 6000 K, and adding the departure enthalpy evaluated using the PR EOS. The reference state for the enthalpy is the triple point of water at 0.01 °C. CO2 properties calculated with EOS2H and REFROP8 (NIST, 2007) are compared in Fig. 2 for temperatures (T) from 200 to 1000 °C and pressures (P) up to 200 MPa. For the properties of the H2O-CO2 mixtures, including the calculation of phases composition from the set of primary variables, the same approach as used in EOS2 is followed, with a few additional options (Battistelli et al., 2017):
supercritical conditions within the family of TOUGH2 reservoir simulators, the existing EOS2 module for H2O-CO2 mixtures (O’Sullivan et al., 1985) has been extended to supercritical gas-like conditions (Saipem, 2016) in order to model the IAPWS (2007) supercritical Region 2 (up to 800 °C) and Region 5 (up to 2000 °C), in addition to the subcritical Regions 1 (liquid water up to 1000 bar), 2 (steam) and 4 (saturation line), which are conventionally modelled by TOUGH2 V.2.0 up to 350 °C. The enhanced module, called EOS2H, has then been coupled (Saipem, 2017) to T2Well (Pan and Oldenburg, 2014), the TOUGH2 extension developed for the transient modeling of fully coupled wellbore-reservoir flow. The resulting code has been verified by comparison with published numerical reservoir solutions and against both experimental and numerical subcritical wellbore flow P&T data. Finally, T2Well-EOS2H has been applied to reproduce the output curve of well IDDP-1 drilled in the Krafla field, Iceland, intercepting a relatively shallow steam-like supercritical geothermal reservoir.
- Under subcritical conditions, the CO2 solubility can be calculated including a Poynting correction in which the molar volume of dissolved CO2 has been used as a regression parameter in order to reproduce the CO2 solubility at 300 bar given by Duan and Sun (2003) from 10 to 250 °C and by Spycher and Pruess (2010) from 250 to 300 °C. The calibrated molar volume (a function of temperature) is then used to account for pressure effects on CO2 solubility in pure water. Comparison between computed and measured CO2 solubility at temperatures up to 300 °C and pressures up to 300 bar showed an acceptable agreement. - The CO2 heat of solution is computed in the original EOS2 module using a regression on data published by Ellis and Golding (1963), which were obtained by using the van ‘t Hoff equation and by deriving the Henry's constant determined with experimental CO2 solubility data. As an option, the heat of solution may also be computed using the same approach but deriving the Henry’s constant used in EOS2 as already made in the EWASG EOS module (Battistelli et al., 1997). - The dissolved CO2 affects the aqueous phase enthalpy. The EOS2 approach is based on the calculation of CO2 enthalpy in the gas phase and on the addition of the CO2 heat of solution in the aqueous
2. EOS2H module and coupling with T2Well EOS2H has been developed by modifying the EOS2 module (O’Sullivan et al., 1985) of the TOUGH2 V.2.0 reservoir simulator, which was already improved as described in Battistelli et al. (2017). EOS2H phase equilibria still use the assumption of low to moderate partial pressures of CO2 (PCO2). As PCO2 is always a primary variable as in EOS2, the partial pressure of water is easily computed and used to i) identify the IAPWS region corresponding to the primary variable set (performed by the new subroutine REGIONS) and ii) to employ the appropriate IAPWS-IF97 correlations for density and enthalpy, and the IAPWS (2008) correlations for viscosity which have been specifically coded for TOUGH2 V.2.0 (Battistelli et al., 2017). In practice, EOS2H uses its own IAPWS-IF97 subroutines and specific phase transition tests which are different than those implemented by Magnusdottir and Finsterle (2015) in EOS1sc. Thus, full verification tests of sub and supercritical gas-like capabilities as well as of coupled wellbore – reservoir flow were needed. Following Magnusdottir and Finsterle (2015), as shown in Fig. 1, in EOS2H the Region 5 correlations are used up to 200 MPa, and Region 2 correlations are used up to 200 MPa above 590 °C. Region 1 correlations
Fig. 1. P&T fields simulated by EOS2H: original IAPWS (2007) in green; extrapolated correlations in cyan. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 2
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 2. Density (a), enthalpy (b) and dynamic viscosity (c) of pure CO2 computed using the PR EOS and the LBC against values computed using REFROP8.
or the recovered condensed steam when simulating the exploitation of a geothermal system. EOS2H has then been coupled (Saipem, 2017) to T2Well (Pan and Oldenburg, 2014), the TOUGH2 extension developed for the transient modeling of fully coupled wellbore-reservoir flow. The T2Well version employed here is the same as that already used with the EWASG EOS module to simulate the coupled flow in geothermal wells (Vasini et al., 2017). This version contains an analytical solution of heat transfer between the wellbore and surrounding rock formations following the approach by Ramey (1962) which includes both the steady-state heat conduction through the components of well completion (casings and cement sheaths) and the radial transient heat conduction within the rock formation. Ramey’s analytical function for heat exchange between wellbore and formation is substituted by the Chiu and Thakur (1991) time function.
phase which is computed neglecting the pressure effects. This approach is acceptable for low partial pressures of CO2, at which its enthalpy in the gas phase can be evaluated with an ideal gas approach. When the partial pressure becomes high, the CO2 gas enthalpy depends also on pressure, and the pressure effect is simply transmitted to the enthalpy of dissolved CO2 if the heat of solution is computed neglecting the pressure. Koschel et al. (2006) showed that the enthalpy of CO2 dissolved in the aqueous phase is almost constant with respect to pressure, indicating that the effect of pressure on the heat of solution is almost compensating for that on gaseous CO2 enthalpy. Thus, in addition of the previous approach a new one can be optionally chosen to calculate the enthalpy of dissolved CO2 as a function of temperature only. - Under sub-critical conditions, phase transitions taking place during TOUGH2 numerical simulations are occasionally affected by numerical instabilities which, among other issues, may be related to the difficulties to interpolate the fluid properties at the interface between grid blocks. The latter problem has been studied by O’Sullivan et al. (2013, 2014) for the EOS modules commonly employed for geothermal reservoir simulation, including EOS2. In order to obtain a more stable interpolation of phase density across a phase transition, they implemented the calculation of the density of non-existing phases. The same approach is now optionally available in EOS2H. In addition, the same option performs also the calculation of enthalpy and viscosity of the non-existing phase, in a way similar to that implemented into iTOUGH2 (Finsterle, 2015).
3. T2well-EOS2H verification 3.1. Verification of reservoir flow against iTOUGH2-EOS1sc The initial verification of EOS2H has been conducted by using the TOUGH2 RFP sample problem (Pruess et al., 1999) involving the injection and production at constant rate of 30 kg/s using a 5-spot well pattern. Taking advantage of the symmetry, the grid is limited to 1/8 of the well pattern and is built parallel to the diagonal connecting injection and production wells. EOS2H has successfully reproduced (Saipem, 2016) (a) the original RFP sample problem results simulated for pure water with EOS1 and EOS2 with reservoir T of 300 °C, (b) the EOS2 results at 300 °C and initial PCO2 of 0.5 and 5 MPa, and (c) the iTOUGH2-EOS1sc results for pure water with initial reservoir
In addition, EOS2H can optionally simulate a traced water component in complete analogy to EOS1 (Pruess et al., 1999). The option can be employed to trace reinjected water, such as the separated brine and/ 3
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 3. P vs T conditions simulated halfway between injection and production wells: comparison between results obtained with iTOUGH2-EOS1sc (symbols and lines) and TOUGH2-EOS2H (lines). Modified after Magnusdottir and Finsterle (2015). IAPWS-IF97 values computed by EOS1sc are overlapped by those computed by EOS2H.
correlation to the other has not shown any numerical problems in the simulated cases: with the maximum time step of 36.5 days, no time step reductions were encountered, nonlinearities were solved with a reasonable number of Newtonian iterations.
conditions of 1100 °C and 90 MPa, 1500 °C and 150 MPa, and 1200 °C and 50 MPa (Magnusdottir and Finsterle, 2015). Comparison of EOS2H results and those simulated with iTOUGH2-EOS1sc are shown in Fig. 3 for the grid element located halfway between the injector and producer. The somewhat unrealistic initial reservoir supercritical conditions (with temperatures above the melting point of commonly encountered rocks) were chosen to replicate the previous iTOUGH2-EOS1sc results, and to show that extrapolated correlations for Regions 2 and 5 do work properly. iTOUGH2-EOS1sc results are correctly replicated by EOS2H. The PT path of selected grid elements of the 3 simulations is reported in Fig. 4 on the PT thermodynamic diagram of pure water, on which the limits of IAPWS-IF97 regions are also drawn. The PT path of production and reinjection wells is shown, as well as that of the element halfway between the two wells. In the T = 1500 °C case, the production well block remains in the extrapolated field of Region 5, while the injection well block goes from the extrapolated field of Region 2 inside Region 3 and then back into the extrapolated field of Region 2. In the T = 1200 °C case, the intermediate element and the injection well block go from the extrapolated field of Region 5 inside Region 2, while the production well block remains in the same Region. In the T = 1100 °C case, the intermediate element and the production well block remain in Region 5, while the injection one enters the extrapolated field of Region 5 to terminate inside Region 2. The simulated cases have involved Regions 2 and 5 and their extrapolated fields above 100 and 50 MPa, respectively, which are the upper limits of the IAPWS-IF97 correlations. The switch from one
3.2. Verification of reservoir flow against AUTOUGH and STAR-HOTH2O Additional simulations have been performed to compare EOS2H with the AUTOUGH results of Croucher and O’Sullivan (2008) who reproduced previous simulations performed with the STAR-HOTH2O simulator (Pritchett, 1994) by Yano and Ishido (1998). A 1D radial model is used with outer radius of 12,039 m, thickness of 100 m, discretized with 40 elements with a radius increment ratio of 1.3. The model features and rock properties listed by Croucher and O’Sullivan (2008) have been used. Two series of simulations were performed with a constant production rate of 15.7 kg/s for 300 h at: A.) reservoir P of 30 MPa and T from 200 to 500 °C; and B.) reservoir T of 450 °C and P from 25 to 35 MPa. The comparison of EOS2H results and those of STAR (black lines) and AUTOUGH (black symbols) are shown in Figs. 5 and 6. Case A at 400 °C (Fig. 5) cannot be modeled with EOS2H as reservoir conditions corresponding to Region 3 are encountered. For the other cases a good reproduction of AUTOUGH results is shown, while the departure from STAR results at early times is due to a difference in near well 4
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 4. P&T thermodynamic diagram of pure water: evolution of reservoir conditions at the injection and production well locations and at the halfway point for the 3 simulated cases.
analysis to the simulated drawdown results. Simulations have been performed starting with 23 MPa initial reservoir pressure and a CO2 partial pressure of 2 MPa, and with initial temperatures of 350, 370, 400, 450 and 500 °C. The partial pressure of water is then equal to 21 MPa, which is slightly lower than the water saturation pressure at 370 °C equal to 21.0434 MPa. Then, the simulation at 370 °C has been performed by assigning two-phase initial conditions with a gas phase saturation of 0.5; the simulation at 350 °C refers to single-liquid conditions; those at 400, 450 and 500 °C refer to supercritical steam-like conditions. In all cases the well is operated at a constant mass extraction rate of 15.7 kg/s for 300 h (1.08 × 106 s). The resulting pressures as a function of time in the well block
permeability introduced by Croucher and O’Sullivan (2008). Fig. 6 also shows a good reproduction of AUTOUGH and STAR results, except for early time P computed with STAR by Yano and Ishido (1998). The comparison confirms that the IAPWS-IF97 correlations have been correctly included into TOUGH2, and that no adverse numerical problems are encountered using EOS2H. 3.3. Code demonstration against a simple drawdown solution The reservoir model used in Section 3.2 above has been employed with different initial reservoir conditions to simulate short production tests of 300 h in order to apply a customary semi-log pressure transient
Fig. 5. Wellbore P vs time: reservoir P 30 MPa. EOS2H vs AUTOUGH (symbols) and STAR (lines) results (modified after Croucher and O’Sullivan, 2008). 5
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 6. P change vs time: reservoir T 450 °C. EOS2H vs AUTOUGH (symbols) and STAR (lines) results (modified after Croucher and O’Sullivan, 2008).
Fig. 7. P = 23 MPa and CO2 partial pressure of 2 MPa. Well block pressure is shown as function of production time on a linear scale.
evaluation of reservoir hydraulic transmissivity (m3) with the following equation:
element are shown in Figs. 7 and 8, on linear and logarithmic scales, respectively. The pressure drawdown at the well block increases with increasing temperature, except for the 370 °C case, as the two-phase conditions lead to a greater compressibility than the single-liquid case at 350 °C. The fluid density decreases with increasing temperatures, and then the drawdown is greater when the extraction rate is held constant. The semi-log plot in Fig. 8 highlights linear trends at late time corresponding to the radial flow dominated period which can be used for the evaluation of permeability-thickness product kh in analogy to the pressure data recorded downhole during a drawdown pressure transient. In Fig. 8 the data at late time regressed with a straight line are shown as function of the logarithm of time. For each case the regression straight line has been determined with the evaluation of the slope m (Pa/cycle) which can be used for the
kh = 2.3026
Wμ 4πρ m
(1)
where: k permeability, m2 h layer thickness, m W mass rate, kg/s μ dynamic viscosity, Pa s ρ density, kg/m3 Tables 1 and 2 list the slope values of the regression straight line, the permeability k evaluated with Eq. 1 using the fluid viscosity and density at i) initial reservoir conditions and ii) in the well block at the 6
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 8. P = 23 MPa and CO2 partial pressure of 2 MPa. Well block pressure is shown as function of logarithm of production time.
changes more at increasing reservoir temperature, and a temperature decline on the order of 10 °C is also simulated as shown in Fig. 9. As viscosity and density of steam and CO2 mixture depend on both pressure and temperature, they change more due to reservoir exploitation. Eq. 1 is strictly valid only for a low compressibility fluid with constant thermo-physical properties. Thus, the difference between the computed permeability value and that used for the numerical simulations is due to the non-negligible variations of supercritical mixture properties during reservoir pressure drawdown.
Table 1 Slope of regression straight line and kh product evaluated with viscosity and density values corresponding to initial reservoir conditions for single-phase cases.
m (Pa/cycle) μ (Pa s) q (m3/s) ρ (kg/m3) kh (m3) k (m2) Error k (%)
350°C
400°C
450°C
500°C
3.42E+05 7.15E-05 2.55E-02 6.17E+02 9.76E-13 9.76E-15 −2.37
7.35E+05 2.67E-05 1.25E-01 1.26E+02 8.32E-13 8.32E-15 −16.81
1.05E+06 2.84E-05 1.59E-01 9.88E+01 7.87E-13 7.87E-15 −21.27
1.40E+06 3.04E-05 1.84E-01 8.55E+01 7.29E-13 7.29E-15 −27.10
3.4. Verification of steady state wellbore flow against T2Well-EWASG and PROFILI simulators T2Well-EOS2H has been compared to T2Well-EWASG by simulating the steady-state wellbore flow in two geothermal wells discharging fluids with high contents of CO2; details are given by Vasini et al. (2018). The simulations are related to:
Table 2 Slope of regression straight line and kh product evaluated with viscosity and density values corresponding to final well block reservoir conditions for singlephase cases.
m (Pa/cycle) μ (Pa s) q (m3/s) ρ (kg/m3) kh /m3) k (m2) Error k (%)
350°C
400°C
450°C
500°C
3.42E+05 7.03E-05 2.58E-02 6.08E+02 9.73E-13 9.73E-15 −2.68
7.35E+05 2.54E-05 1.60E-01 9.81E+01 1.02E-12 1.02E-14 +1.58
1.05E+06 2.73E-05 2.18E-01 7.22E+01 1.03E-12 1.03E-14 +3.40
1.40E+06 2.93E-05 2.72E-01 5.78E+01 1.04E-12 1.04E-14 +3.79
a) Well KD-13, Kizildere field, Turkey, discharging a low NaCl brine (1000 ppm) with a high content of dissolved CO2 (20,000 ppm); data published by James (1975); b) Well W2 from an unspecified field managed by ENEL, discharging a brine with remarkable NaCl (9600 ppm) and CO2 (30,000 ppm) contents; data published by Barelli et al. (1982). The results of the comparison are shown in Fig. 10, where field measurements are presented together with T2Well-EWASG and T2WellEOS2H results. Even though slight differences are present, which are most likely related to a different wellbore gridding and to the need to adjust the CO2 concentration due to the absence of NaCl effects, the results show that T2Well-EOS2H is able to reproduce both the experimental steady-state wellbore flow data and the corresponding T2WellEWASG simulation results of Vasini et al. (2018). T2Well-EOS2H is also compared to a steady-state wellbore flow simulation performed with PROFILI (Saipem, 2017) based on conditions of Well IDDP-1, Krafla field, Iceland (Einarsson et al., 2015). The numerical wellbore flow simulator PROFILI (Battistelli, 2010) simulates the steady-state flow under production or injection of H2ONaCl-CO2 fluid mixtures under single-phase liquid, two-phase and single-phase gas conditions (including supercritical gas-like conditions)
end of 300 h of discharge, respectively. While the permeability for the 350 °C single-phase liquid is evaluated with a comparable error for the different fluid parameters (−2.37 and −2.68 %, respectively), for the cases with supercritical temperature, the error on permeability is much higher, between −16 and −27 % when using the initial reservoir conditions, with respect to that obtained using the final well block fluid parameters, which is between 1.58 and 3.79 %. For the single-phase liquid case at 350 °C, no phase transitions are experienced and then both fluid viscosity and density remain almost constant as the pressure transient occurs at almost isothermal conditions and liquid fluid properties do not change considerably with pressure. For supercritical steam-like conditions, the well block pressure 7
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 9. Well block P&T path for the 5 simulated cases: reservoir pressure of 23 MPa and CO2 partial pressure of 2 MPa.
4. Modeling of coupled wellbore-reservoir flow at well sector scale
in vertical or directional wellbores with variable diameters and hole surface roughness. The numerical solution of mass, energy and momentum balance equations is based on the work by Barelli et al. (1982). For two-phase conditions, pressure changes due to friction, acceleration and gravity can be evaluated using a homogeneous approach based on the DIF-3 correlations by CISE (Lombardi and Ceresa, 1978; Bonfanti et al., 1979) or the CESNEF-2 correlations (Lombardi and Carsana, 1992). Well IDDP-1 is completed with a 9-5/8” casing down to 1950 m and with a 9-5/8” slotted liner down to 2070 m. The wellbore is here assumed to be vertical even though the production tests were performed using the final wellbore path resulting from two side-tracks. The PROFILI simulation is performed from well-head to bottom-hole. The rate is 7 kg/s with wellhead P (WHP) of 14.2 MPa, wellhead T (WHT) of 440 °C, and 1000 ppm CO2. Heat transfer is computed following Ramey (1962) with an assumed average global heat exchange coefficient U of 20 W/(m2 °C). The T2Well-EOS2H simulation is performed by fixing flowing bottom conditions equal to those computed by PROFILI and extracting 7 kg/s from the well-head until almost steady-state conditions are reached after 30 days. Fig. 11 shows good agreement between flowing P&T computed by PROFILI and T2Well-EOS2H.
A preliminary coupled wellbore-reservoir flow simulation was previously presented by Battistelli et al. (2018) based on the results of production tests performed at well IDDP-1, Krafla, Iceland (here referred as 2018 results). The simulations did not intend to accurately replicate the IDDP-1 production test results as testing of the well was complex and performed in 5 stages of quite different duration from March 2010 to June 2012. The production testing assemblage was changed according to encountered operating conditions with some long shut-in periods between production stages. In addition, Einarsson et al. (2015) provide the WHP, WHT and wellhead enthalpy (WHE) history as functions of time but not the related discharge rate history. Thus, because of the complexity of IDDP-1 testing operations and the lack of some important data, the simulations did not intend to accurately reproduce the observed well behavior, but just consider IDDP-1 published results as a realistic framework to be used for the modeling of coupled wellbore-reservoir flow. First, the well output curve (OC) reconstructed by Einarsson et al. (2015) was reproduced using PROFILI by adjusting the static reservoir P&T and using a parabolic P draw-down as a function of mass rate with
Fig. 10. Simulated P&T vs depth: (a) well KD13; (b) well W2. T2Well-EOS2H vs T2Well-EWASG. 8
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 11. Simulated P&T vs depth for well IDDP-1: T2Well-EOS2H vs PROFILI comparison.
while both WHT and WHE are slightly higher due to the higher overall heat exchange coefficient of well completion with respect to the average value used previously. The OC is then reproduced with T2Well-EOS2H by simulating in full transient conditions well discharge at increasing rates from 5 to 50 kg/s with constant production steps lasting 6 days each as in the new PROFILI steady state simulation. The reservoir is described using the same 1D radial model of 1 km radius and 100 m thickness (from 2000 to 2100 m depth) previously employed. The results in terms of WHP and WHT for 2018 and 2019-A simulations are shown in Fig. 13; they are obtained by trial and error using permeabilities of 8.45 × 10−14 m2 and 8.70 × 10−14 m2, respectively. Porosity has been set at 0.10 and no skin effects have been included. WHP and WHT are plotted vs discharge time, while the steady-state values computed using PROFILI are shown as circles at the end of each production step. WHP was previously (Battistelli et al., 2018) slightly underestimated by T2Well-EOS2H for rates of 15–30 kg/s and overestimated at the highest rate of 50 kg/s, when the fluid velocity simulated by PROFILI reached 209 m/s at the well-head, being only 51 m/s at bottom-hole. At simulated wellhead conditions the fluid velocity was still much lower than the speed of sound of pure water which is about 646 m/s. The friction gradient increased from 2.50 × 103 Pa/m at bottom-hole to 7.05 × 103 Pa/m at the wellhead. At 45 kg/s the fluid velocity increased from bottom to well-head from 40.6 m/s to only 67.0 m/s. Thus, the increment from 45 to 50 kg/s produced an important increment of fluid velocity which was nevertheless still well below chocked flow conditions. The new 2019-A simulation shows similar trends for WHP and WHT, with WHP lower than PROFILI at rates lower than 40 kg/s, and higher at rates above 40 kg/s. The increment of production step duration from 3 to 6 days induced a non-negligible inflow of supercritical fluid from the lateral outer boundary at the highest rates. In order to have an almost infinite acting reservoir, a new grid was built with the lateral boundary set at 5 km, and with the same discretization as the previous grid in the first 1000 m. The BHP and BHT computed with T2Well-EOS2H (called 2019B) are also shown in Fig. 13. They were obtained by increasing the reservoir permeability to 9.30 × 10−14 m2. Apart from that, the results of the simulation 2019-A are very close to those of simulation 2019-B. The latter results are used in the subsequent comparison with the PROFILI 2019 results.
viscous and turbulent coefficients of 0.5 bar/(kg/s) and 0.0175 bar/ (kg/s)2, respectively. The global heat exchange coefficient of well completion, U, was assumed constant at 10 W/(m2 °C) to compute the heat exchange with surrounding rock formations following Ramey (1962). The heat conductivity and diffusivity of the surrounding rock formations were assumed constant with depth and set to 2.25 W/(m °C) and 9.3E-7 m2/s, respectively. The IDDP-1 production casing is a Hydrill 563, 9-5/8″ 53.5 lb/ft pipe, while the slotted liner is a 9-5/8″ 47 lb/ft pipe with Buttress connections (Pálsson et al., 2014). The roughness of casing and liner was assumed to be 1.0E-5 m and 1.0E4 m, respectively. The output curve was computed at increasing rates from 5 to 50 kg/s, assuming a constant average production time of 30 days for each rate. While the bottom hole P (BHP) changes with rate, the bottom hole T (BHT) was held constant for each production step, and a PCO2 of 7.2 kPa corresponding to 1000 ppm was assumed. Initial static P&T values were calibrated at 15.5 MPa and 510 °C. The OC was then reproduced with T2Well-EOS2H by simulating well discharge at increasing rates from 5 to 50 kg/s with constant production steps lasting only 3 days each, instead of 30 days. The reservoir was described using a 1D radial model of 1 km radius and 100 m thickness (from 2000 to 2100 m depth) with homogeneous, isotropic rock properties, initial P&T and fluid composition equal to that calibrated using the PROFILI simulation. The radial grid includes the wellbore block plus 100 radial elements with the lateral boundary at 1 km kept at its initial conditions. All the parameters related to the heat exchange were the same as those used for PROFILI simulations. For the new PROFILI simulation described below (referred to as 2019 results), the global heat exchange coefficient of well completion, U, is computed following Willhite (1967) and using the final well design given by Pálsson et al. (2014) at 10.25 W/(m2 °C) down to 254 m, 12.19 W/(m2 °C) down to 785 m, and 26.86 W/(m2 °C) down to 1900 m. The output curve is computed at increasing rates from 5 to 50 kg/s, with 10 steps of 6 days each and the production time for each step is varied accordingly. The same parameters describing the downhole deliverability curve already calibrated in Battistelli et al. (2018) were used. Fig. 12 shows the 2018 and 2019 results of PROFILI simulations with the OC and WHT on the left, and the WHE on the right. Changes made in the calculation of heat exchange had minor effects on the OC which is almost unchanged with respect to the previous simulation, 9
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 12. IDDP-1 output curve: mass rate and WHT vs WHP (a) and mass rate and WHE vs WHP (b) simulated using PROFILI (modified after Einarsson et al., 2015). “Calculated mass Flow” and “Characteristic Curve” are the experimental rates and interpreted output curve, respectively, of Einarsson et al. (2015).
high rates. In fact, water specific heat at 15 MPa and T from 410 to 510 °C changes from 3.91 to 2.85 kJ/(kg °C), while it is substantially higher for subcritical liquid water as it ranges between 4.15 and 5.47 kJ/(kg °C) from 50 to 300 °C. Geothermal wells discharging supercritical fluids are likely to experience high heat losses in the relatively cold upper well section which turns into a relevant T reduction of discharged fluid due to its low heat capacity. Among the reasons for the WHT differences between PROFILI and T2Well we can first mention that while PROFILI uses a steady-state approach, T2Well computes a transient solution even when the heat transfer is computed with the analytical method by Ramey (1962). Moreover, the PROFILI simulations have been performed by assuming a constant BHT, as the simulator just computes the wellbore flow. On the other hand, T2Well-EOS2H computes the BHT as a result of reservoir flow by solving both mass and heat balance equations. Fig. 16 shows
A direct comparison between the 2019 output curve simulated by PROFILI and T2Well-EOS2H (2019-B) is shown in Fig. 14, where both mass rate and WHT are plotted versus WHP. Fig. 15 shows WHP and WHT versus mass rate. The T2Well-EOS2H WHP is lower and higher than PROFILI above and below 40 kg/s, respectively. The maximum discrepancy is obtained at the highest rate of 50 kg/s where the WHP computed by T2Well-EOS2H is significantly higher. Fig. 15 shows how the WHT increases considerably with rate up to about 15 kg/s mostly because the higher rates imply smaller heat losses per unit mass. Then the WHT almost levels out up to about 40 kg/s. It finally declines at the highest rates. Differences between the two sets of results is within ± 5 °C, apart for the 45 kg/s results. Temperature changes for wells discharging supercritical water seem to be relatively sensitive to enthalpy variations, as those related to the heat exchange with surrounding formation and kinetic energy losses at
Fig. 13. IDDP-1 output curve: WHP and WHT vs time simulated using T2Well-EOS2H. 2018 (a), 2019-A (b) and 2019-B (c) simulation results compared with corresponding steady-state PROFILI results plotted at the end of each step. 10
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 14. IDDP-1 output curve: Mass rate and WHT vs WHP simulated using T2Well-EOS2H, compared with values computed by PROFILI and the experimental ‘Characteristic Curve’ by Einarsson et al. (2015).
Fig. 15. IDDP-1 output curve: WHP and WHT vs mass rate simulated using T2Well-EOS2H, compared with values computed by PROFILI.
could be a better assumption for the PROFILI wellbore flow steady-state simulations. Then, the OC was simulated again using PROFILI but assuming a constant BHE equal to the initial reservoir value. The computed BHT as function of rate is shown in Fig. 18 which suggests that with the assumption of a constant BHE the BHT are lower than those actually computed by T2Well-EOS2H. Fig. 19 shows the effect of the constant BHE assumption on the PROFILI output curve: the computed WHT becomes much lower than the previous 2019 results, while the effects on the WHP are significant only at the highest rate of 50 kg/s. In practice, because of the strong dependency of supercritical steam enthalpy on both P&T, the common assumption of assigning the BH conditions in wellbore flow simulation
the BHP and BHT as function of rate for both PROFILI and T2Well. The BHT computed by T2Well declines at increasing rates, becoming 57.7 °C lower at the highest rate of 50 kg/s. As can be expected, the enthalpy at bottom hole has also a different trend for the two simulations. Starting with the same value of 3330 kJ/ kg at initial reservoir conditions, PROFILI computes a final BHE of 3414 kJ/kg greater than the 3359 kJ/kg of T2Well, because of the PROFILI assumption of constant BHT. Fig. 17 shows that the WHE increases at increasing rates up to the maximum simulated rate of 50 kg/s. The increment is higher up to 10−15 kg/s. Then the WHE increment is reduced as the rate increases. The BHE history of T2Well suggests checking if a constant BHE
Fig. 16. IDDP-1 downhole deliverability curve: BHP and BHT vs mass rate simulated using T2Well-EOS2H, compared with values used by PROFILI. 11
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 17. IDDP-1 output curve: WHP and WHE vs mass rate simulated using T2Well-EOS2H, compared with values computed by PROFILI.
Fig. 18. IDDP-1 downhole deliverability curve: BHP and BHT vs mass rate simulated using T2Well-EOS2H, (OC 2019-B) compared with BHT values computed by PROFILI when assuming constant BHE.
EOS2H is capable of accurately modeling coupled wellbore-reservoir flow processes under supercritical steam-like conditions.
for wells producing from liquid dominated reservoirs do not simply apply when modeling the flow of supercritical fluids. The coupled wellbore-reservoir flow modeling seems to be able to supply more reliable results than a wellbore flow simulator. While the simulated temperatures might be considerably affected by the accuracy of heat balance calculations and related assumptions, computed flowing pressures seem to be less affected. These preliminary simulations of coupled wellbore-reservoir flow at a well sector scale show promising results, suggesting that T2Well-
5. Conclusions The existing EOS2 module of the TOUGH2 V.2.0 reservoir simulator for H2O-CO2 mixtures has been improved by including the capability to model the thermodynamic conditions of supercritical steam-like reservoirs using the IAPWS-IF97 correlations within Regions 2 and 5. The
Fig. 19. IDDP-1 output curve: Mass rate and WHT vs WHP simulated using T2Well-EOS2H and the experimental output curve (Einarsson et al., 2015), compared with values computed by PROFILI using both constant BHT and BHE assumptions. 12
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Validation, Writing - original draft. Stefan Finsterle: Methodology, Writing - review & editing. Marica Marcolini: Writing - review & editing, Project administration. Lehua Pan: Methodology, Writing review & editing.
calculation of CO2 properties at high T and P has been improved by including the use of PR EOS for density and enthalpy departure as well as the LBC method for the dynamic viscosity. The improved module, called EOS2H, has been linked to T2Well for the modeling of coupled wellbore-reservoir flow under fully transient conditions in geothermal reservoirs under both sub-critical and supercritical steam-like conditions (Regions 2 and 5 of the IAPWS-IF97 correlations). T2Well-EOS2H accuracy and numerical performances have been successfully verified by comparison with:
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments
- TOUGH2-EOS2 V.2.0 reservoir simulations at sub-critical conditions; - iTOUGH2-EOS1sc reservoir simulations at supercritical conditions without CO2; - STAR-HOTH2O and AUTOUGH reservoir simulations at supercritical conditions; - T2Well-EWASG and PROFILI wellbore flow simulations.
This work was supported by Saipem SpA under the R&D project "Modeling of flow in geothermal wells", Job P5000Q-06. The management of Saipem SpA is acknowledged for the permission to publish this paper. References
The verified code has then been used to model the output curve of a well producing from a supercritical steam-like geothermal reservoir. Wellbore and reservoir characteristics are based on those of well IDDP1 drilled in the Krafla field, Iceland. The output curve of well IDDP-1 inferred from discharge tests has been replicated using the PROFILI steady-state numerical wellbore simulator with capabilities to simulate steam-like supercritical flow conditions, by calibrating the coefficients of the downhole deliverability curve. Then the output curve has been reproduced with T2Well-EOS2H using a simple model with a 1D vertical grid for the wellbore and a 1D radial grid for the supercritical reservoir. The numerical model calibration has been limited to a simple trial and error approach by modifying reservoir permeability. A more rigorous approach was deemed to be not necessary because several of the wellbore and reservoir features used were not available and then were just assumed based on educated guesses. T2Well-EOS2H simulated output curve shows differences with PROFILI results which are likely to be due to:
Barelli, A., Corsi, R., Del Pizzo, G., Scali, C., 1982. A two-phase flow model for geothermal wells in the presence of non-condensable gas. Geothermics 11 (3), 175–191. Battistelli, A., Calore, C., Pruess, K., 1997. The simulator TOUGH2/EWASG for modelling geothermal reservoirs with brines and a non-condensible gas. Geothermics 26 (4), 437–464. Battistelli, A., 2010. PROFILI Wellbore Flow Simulator: From Geothermal Wells to GHG and Acid Gas Injection Wells. GeoThermExpo 2010, Ferrara September 21-22, 2010. Battistelli, A., Swenson, D., Alcott, A., 2017. Improved PetraSim-TOUGH2 Capabilities for the Simulation of Geothermal Reservoirs. 42nd Stanford Geoth. Work, Feb. 13-15, 2017, Stanford, CA. Battistelli, A., Finsterle, S., Pan, L., 2018. Modeling of coupled wellbore-reservoir flow for supercritical geothermal systems in steam-like conditions. In: TOUGH Symposium 2018. Berkeley, CA. Bonfanti, F., Ceresa, I., Lombardi, C., 1979. Two-phase pressure drops in the low flowrate region. Energia Nucleare 26 (10), 481–492. Brikowski, T.H., 2001. Modeling supercritical systems with TOUGH2: preliminary results using the EOS1SC equation of state module. In: Proc. 26th Work. Geoth. Res. Eng., Stanford University. Stanford, CA. Chase Jr, M.W., 1998. NIST-JANAF themochemical tables. J. Phys. Chem. Ref. Data, Monograph 9, Fourth Edition. pp. 1–1951. Chiu, K., Thakur, S.C., 1991. Modeling of wellbore heat losses in directional wells under changing injection conditions. In: SPE 22870. Presented at SPE Annual Technical Conference and Exhibition. October, 1991. Dallas, Texas. pp. 517–528. Cox, B.L., Pruess, K., 1990. Numerical experiments on convective heat transfer in watersaturated porous media at near-critical conditions. Transp. Porous Media 5 (3), 299–323. Croucher, A.E., O’Sullivan, M.J., 2008. Application of the computer code TOUGH2 to the simulation of supercritical conditions in geothermal systems. Geothermics 37, 622–634. Duan, Z.H., Sun, R., 2003. An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chem. Geol. 193, 257–271. Einarsson, K.K., Sveinsson, E., Ingason, K., Kristjánsson, V., Hólmgeirsson, S., 2015. Discharge testing of magma Well IDDP-1. In: Proc. WGC 2015. Melbourne, Australia.. Ellis, A.J., Golding, R.M., 1963. The solubility of carbon dioxide above 100°C in water and in sodium chloride solutions. Am. J. Sci. 261, 47–60. Finsterle, S., Commer, M., Edmiston, J., Jung, Y., Kowalsky, M.B., Pau, G.S.H., Wainwright, H., Zhang, Y., 2017. iTOUGH2: a multiphysics simulation-optimization framework for analyzing subsurface systems. Comput. Geosci. 108, 8–20. Finsterle, S., 2015. Enhancements to the TOUGH2 Simulator Implemented in iTOUGH2, Report LBNL-7016E. Lawrence Berkeley National Laboratory, Berkeley, CA. IAPWS, 2007. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam. Switzerland. IAPWS, 2008. Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance. Berlin, Germany. James, R., 1975. Gas content of a hot-water reservoir estimated from downhole pressure and temperature measurements. In: Berkeley, CA. Second United Nations Symposium Vol. 3. pp. 1689–1691. Kissling, W.M., 1995. Extending MULKOM to super-critical temperatures and pressures. In: Proc. WGC 1995. Firenze, Italy. pp. 1687–1690. Kissling, W.M., 2005. Transport of three-phase hyper-saline brines in porous media: theory and code implementation. Transp. Porous Media 61, 25–44. Koschel, D., Coxam, J.-Y., Rodier, L., Majer, V., 2006. Enthalpy and solubility data of CO2 in water and NaCl(aq) at conditions of interest for geological sequestration. Fluid Phase Equilib. 247, 107–120. Lohrentz, J., Bray, B.G., Clark, C.R., 1964. Calculating viscosities of reservoir fluids from their compositions. SPE paper 915. J. Pet. Technol. Altern. Fuels 1171–1176. Lombardi, C., Ceresa, I., 1978. A generalized pressure drop correlation in two-phase flow. Energia Nucleare 25 (4), 181–198. Lombardi, C., Carsana, C.G., 1992. A dimensionless pressure drop correlation for two-
i) the different numerical solutions implemented, with steady-state wellbore flow for PROFILI, and full transient wellbore-reservoir flow for T2Well-EOS2H, and ii) the use of a downhole deliverability curve for PROFILI with the assumption of fixing bottom hole flowing conditions as either constant BHT or constant BHE, while variable BHT and BHE are actually computed by T2Well-EOS2H. The simulations of wellbore flow in the presence of supercritical fluids suggest that an accurate modeling of heat transfer toward the surrounding rock formation might be more important than for conventional sub-critical geothermal wells. In fact, geothermal wells discharging supercritical fluids are likely to experience large heat losses in the relatively cold upper well section which turns into a significant temperature reduction of the discharged fluid due to its low heat capacity. This aspect requires additional efforts in order to improve the present modeling capabilities, even though using a full numerical approach for the modeling of radial heat conduction would probably be already able to improve the reliability of the results. The preliminary simulations of coupled wellbore-reservoir flow presented above show promising results, suggesting that T2Well-EOS2H is capable of accurately modeling the coupled wellbore-reservoir flow under supercritical steam-like conditions as those encountered near magmatic heat sources and within the deep roots of already exploited geothermal fields. CRediT authorship contribution statement Alfredo Battistelli: Conceptualization, Methodology, Software, 13
Geothermics 86 (2020) 101793
A. Battistelli, et al. phase mixtures flowing up in vertical duct covering wide parameter ranger. Heat Pipe Sci. Technol. Int. J. 10 (1-2), 125–141. Magnusdottir, L., Finsterle, S., 2015. An iTOUGH2 equation-of-state module for modeling supercritical conditions in geothermal reservoirs. Geothermics 57, 8–17. McKibbin, R., McNabb, A., 1999. Deep hydrothermal systems: mathematical modeling of hot dense brines containing noncondensible gases. J. Porous Media 2 (1), 107–126. Montegrossi, G., Pasqua, C., Battistelli, A., Mwawongo, G., Ofwona, C., 2015. 3D natural State model of the menengai geothermal system, Kenya. In: Proc. WGC 2015. Melbourne, Australia, 19–25 April 2015. Morin, A., Lövfall, B.T., Meese, E., 2016. Simulation of Supercritical Water Flow in the Venelle 2 Well. (Accessed 24 June 2019). http://www.descramble-h2020.eu/ images/download/publications/GeoThArt.pdf. NIST, 2007. Standard Reference Database 23. NIST Thermophysical Properties of Hydrocarbon Mixtures – Program REFPROP – Version 8.0.Standard Reference Database 23. NIST Thermophysical Properties of Hydrocarbon Mixtures – Program REFPROP – Version 8.0. O’Sullivan, J., Croucher, A., Yeh, A., O’Sullivan, M.J., 2013. Improved convergence for air-water and CO2-water TOUGH2 simulations. In: Proceedings 35th New Zealand Geoth. Work. Rotorua, New Zealand. O’Sullivan, J., Croucher, A., Yeh, A., O’Sullivan, M.J., 2014. Further improvements in the convergence of TOUGH2 simulations. In: Oñate, E., Oliver, J., Huerta, A. (Eds.), 11th World Congress on Computational Mechanics (WCCM XI). O’Sullivan, J., Kipyego, E., Croucher, A., Ofwona, C., O’Sullivan, M.J., 2015. A supercritical model of the menengai geothermal system. In: Proc. WGC 2015. Melbourne, Australia, 19–25 April 2015. O’Sullivan, J., O’Sullivan, M.J., Croucher, A., 2016. Improvements to the AUTOUGH2 supercritical simulator with extension to the air-water equation-of-State. GRC Transactions 40, 921–930. O’Sullivan, M.J., Bodvarsson, G.S., Pruess, K., Blakeley, M.R., 1985. Fluid and heat flow in
gas-rich geothermal reservoirs. SPE J. 25, 215–226. Pálsson, B., Hólmgeirsson, S., Guðmundsson, Á., Bóasson, H.Á., Ingason, K., Sverrison, H., Thórhallsson, S., 2014. Drilling of the well IDDP-1. Geothermics 49, 23–30. Pan, L., Oldenburg, C.M., 2014. T2Well - An integrated wellbore–reservoir simulator. Comp. Geosci. 65, 46–55. Peng, D.-Y., Robinson, D.B., 1976. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 15, 59–64. Pritchett, J.W., 1994. HOTH2O A Description of H2O Properties to one Kilobar and 800°C for Use With the STAR Geothermal Reservoir Simulator. SSS-TR-94-14730, S-Cubed, La Jolla. Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User’s Guide, Version 2.0, Report LBNL-43134. Lawrence Berkeley National Laboratory, Berkeley, CA. Ramey Jr, H.J., 1962. Wellbore heat transmission. In: 36Th Ann. Fall Meet. of SPE, Oct. 811, 1961. Dallas. Saipem, 2016. Numerical Modeling of Hydrothermal Systems with Supercritical Steam Conditions. Module EOS2H, Report SPC_94060, Job P5000Q-06_03 (unpublished). . Saipem, 2017. Simulation of Production Tests With TOUGH2-T2Well, Report SPC_94070, P5000Q-06_03 (unpublished). . Spycher, N., Pruess, K., 2010. A phase-partitioning model for CO2–brine mixtures at elevated temperatures and pressures: application to CO2-Enhanced geothermal systems. Transp. Porous Media 82, 173–196. Vasini, E.M., Battistelli, A., Berry, P., Bonduà, S., Bortolotti, V., Cormio, C., Pan, L., 2018. Interpretation of production tests in geothermal wells with T2Well-EWASG. Geothermics 73, 158–167. Willhite, G.P., 1967. Over-all heat transfer coefficients in steam and hot water injection wells. J. Pet. Technol. Altern. Fuels 607–615 Paper SPE 1449. Yano, Y., Ishido, T., 1998. Numerical investigation of production behavior of deep geothermal reservoirs at supercritical conditions. Geothermics 27 (5/6), 705–721.
14
Contents lists available at ScienceDirect
Geothermics journal homepage: www.elsevier.com/locate/geothermics
Modeling of coupled wellbore-reservoir flow in steam-like supercritical geothermal systems
T
Alfredo Battistellia,*, Stefan Finsterleb, Marica Marcolinic, Lehua Pand a
Independent Consultant, 60033, Chiaravalle, AN, Italy Finsterle GeoConsulting, Kensington, CA, 94708, USA c PRCFA Dept., Saipem SpA, Via Toniolo 1, 61032 Fano, PU, Italy d Energy Geosciences Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA, 94720, USA b
A R T I C LE I N FO
A B S T R A C T
Keywords: Coupled wellbore-reservoir flow Geothermal reservoirs iTOUGH2 Supercritical conditions TOUGH2 T2Well
The capabilities to efficiently model the transport of mass and heat in geothermal reservoirs at supercritical conditions are necessary to support the efforts to explore and utilize supercritical hydrothermal systems, as well as to extend the utilization of existing geothermal fields at greater depths, where supercritical temperatures might be encountered. Similar capabilities are needed to model the flow in wells producing supercritical fluids, possibly in a coupled way with the reservoir flow. As a step in this direction, the new EOS2H module for subcritical and steam-like super-critical H2O-CO2 mixtures has been linked to T2Well, the extension of the TOUGH2 numerical reservoir simulator developed for the transient modeling of fully coupled wellbore-reservoir flow. The resulting code is verified against simulations performed with other supercritical reservoir simulators and steadystate subcritical and supercritical wellbore flow simulations. It is then used to perform coupled wellbore-reservoir flow simulations at the well-sector scale relevant to the conditions found in well IDDP-1, Krafla, Iceland.
1. Introduction Supercritical fluids exist near magmatic heat sources, and supercritical conditions are often encountered beneath presently exploited high-temperature geothermal reservoirs. These supercritical ‘roots’ are often neglected when building 3D numerical models of exploited geothermal systems. One example is given by the 3D numerical model of the Menengai geothermal systems (Rift Valley, Kenya) built using TOUGH2-EOS2 (Pruess et al., 1999), which was artificially cut in the up-flow area at about 2000 m below ground level to avoid supercritical temperatures beneath the presently exploited two-phase steam dominated reservoir (Montegrossi et al., 2015). On the other hand, O’Sullivan et al. (2015) presented natural-state simulations at Menengai using a deeper model made possible by the supercritical version of the AUTOUGH simulator (Croucher and O’Sullivan, 2008). Supercritical Equations of State (EOS) for the TOUGH2 and iTOUGH2 (Finsterle et al., 2017) numerical reservoir simulators have been developed in the past considering pure water by Cox and Pruess (1990); Kissling (1995); Brikowski (2001), and more recently by Croucher and O’Sullivan (2008) and Magnusdottir and Finsterle (2015). Developments for mixtures have been presented by Kissling (2005) for
binary mixtures of H2O-NaCl up to 620 °C, and by McKibbin and McNabb (1999) for ternary mixtures of H2O-NaCl-CO2. O’Sullivan et al. (2016) developed an AUTOUGH-based simulator for water and air mixtures in order to incorporate shallow unsaturated zones. Modeling of coupled wellbore-reservoir flow in supercritical conditions is also required to interpret well testing operations, and in general to model wellbore flow effects on the performances of wells discharging from supercritical steam-like reservoirs. Morin et al. (2016) presented steady-state and transient simulations of wellbore flow inspired by the conditions expected for well Venelle 2 in the Larderello field, Italy, which was deepened within the DESCRAMBLED EU sponsored project aimed at experimenting with novel drilling technology and at investigating supercritical conditions in continental Europe. They used the commercial multiphase pipe flow simulator LedaFlow, employed also to simulate the flow in oil and gas production wells. Simulations performed included the heat exchange with surrounding rocks, assuming fixed bottomhole pressure and temperature (P&T) conditions. By simulating start-up and well shut-in operations, they observed that pressures, temperatures and flow rates are considerably affected by the thermal state of the well. As a first step in modeling coupled wellbore-reservoir flow in
⁎
Corresponding author at: AMBEN Dept., Saipem SpA, Via Toniolo 1, 61032 Fano, PU, Italy. E-mail addresses: [email protected] (A. Battistelli), stefan@finsterle-geoconsulting.com (S. Finsterle), [email protected] (M. Marcolini), [email protected] (L. Pan). https://doi.org/10.1016/j.geothermics.2019.101793 Received 25 June 2019; Received in revised form 22 November 2019; Accepted 18 December 2019 0375-6505/ © 2019 Published by Elsevier Ltd.
Geothermics 86 (2020) 101793
A. Battistelli, et al.
can also be used up to 140 MPa and 365 °C with minor errors (Battistelli et al., 2017). Density of pure CO2 is computed using the Peng and Robinson (1976) cubic EOS (PR EOS), while its dynamic viscosity is estimated using the so-called LBC method of Lohrentz et al. (1964). Pure CO2 enthalpy is computed by evaluating the ideal gas enthalpy with the Shomate equation (Chase, 1998), which is valid between 298 and 6000 K, and adding the departure enthalpy evaluated using the PR EOS. The reference state for the enthalpy is the triple point of water at 0.01 °C. CO2 properties calculated with EOS2H and REFROP8 (NIST, 2007) are compared in Fig. 2 for temperatures (T) from 200 to 1000 °C and pressures (P) up to 200 MPa. For the properties of the H2O-CO2 mixtures, including the calculation of phases composition from the set of primary variables, the same approach as used in EOS2 is followed, with a few additional options (Battistelli et al., 2017):
supercritical conditions within the family of TOUGH2 reservoir simulators, the existing EOS2 module for H2O-CO2 mixtures (O’Sullivan et al., 1985) has been extended to supercritical gas-like conditions (Saipem, 2016) in order to model the IAPWS (2007) supercritical Region 2 (up to 800 °C) and Region 5 (up to 2000 °C), in addition to the subcritical Regions 1 (liquid water up to 1000 bar), 2 (steam) and 4 (saturation line), which are conventionally modelled by TOUGH2 V.2.0 up to 350 °C. The enhanced module, called EOS2H, has then been coupled (Saipem, 2017) to T2Well (Pan and Oldenburg, 2014), the TOUGH2 extension developed for the transient modeling of fully coupled wellbore-reservoir flow. The resulting code has been verified by comparison with published numerical reservoir solutions and against both experimental and numerical subcritical wellbore flow P&T data. Finally, T2Well-EOS2H has been applied to reproduce the output curve of well IDDP-1 drilled in the Krafla field, Iceland, intercepting a relatively shallow steam-like supercritical geothermal reservoir.
- Under subcritical conditions, the CO2 solubility can be calculated including a Poynting correction in which the molar volume of dissolved CO2 has been used as a regression parameter in order to reproduce the CO2 solubility at 300 bar given by Duan and Sun (2003) from 10 to 250 °C and by Spycher and Pruess (2010) from 250 to 300 °C. The calibrated molar volume (a function of temperature) is then used to account for pressure effects on CO2 solubility in pure water. Comparison between computed and measured CO2 solubility at temperatures up to 300 °C and pressures up to 300 bar showed an acceptable agreement. - The CO2 heat of solution is computed in the original EOS2 module using a regression on data published by Ellis and Golding (1963), which were obtained by using the van ‘t Hoff equation and by deriving the Henry's constant determined with experimental CO2 solubility data. As an option, the heat of solution may also be computed using the same approach but deriving the Henry’s constant used in EOS2 as already made in the EWASG EOS module (Battistelli et al., 1997). - The dissolved CO2 affects the aqueous phase enthalpy. The EOS2 approach is based on the calculation of CO2 enthalpy in the gas phase and on the addition of the CO2 heat of solution in the aqueous
2. EOS2H module and coupling with T2Well EOS2H has been developed by modifying the EOS2 module (O’Sullivan et al., 1985) of the TOUGH2 V.2.0 reservoir simulator, which was already improved as described in Battistelli et al. (2017). EOS2H phase equilibria still use the assumption of low to moderate partial pressures of CO2 (PCO2). As PCO2 is always a primary variable as in EOS2, the partial pressure of water is easily computed and used to i) identify the IAPWS region corresponding to the primary variable set (performed by the new subroutine REGIONS) and ii) to employ the appropriate IAPWS-IF97 correlations for density and enthalpy, and the IAPWS (2008) correlations for viscosity which have been specifically coded for TOUGH2 V.2.0 (Battistelli et al., 2017). In practice, EOS2H uses its own IAPWS-IF97 subroutines and specific phase transition tests which are different than those implemented by Magnusdottir and Finsterle (2015) in EOS1sc. Thus, full verification tests of sub and supercritical gas-like capabilities as well as of coupled wellbore – reservoir flow were needed. Following Magnusdottir and Finsterle (2015), as shown in Fig. 1, in EOS2H the Region 5 correlations are used up to 200 MPa, and Region 2 correlations are used up to 200 MPa above 590 °C. Region 1 correlations
Fig. 1. P&T fields simulated by EOS2H: original IAPWS (2007) in green; extrapolated correlations in cyan. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 2
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 2. Density (a), enthalpy (b) and dynamic viscosity (c) of pure CO2 computed using the PR EOS and the LBC against values computed using REFROP8.
or the recovered condensed steam when simulating the exploitation of a geothermal system. EOS2H has then been coupled (Saipem, 2017) to T2Well (Pan and Oldenburg, 2014), the TOUGH2 extension developed for the transient modeling of fully coupled wellbore-reservoir flow. The T2Well version employed here is the same as that already used with the EWASG EOS module to simulate the coupled flow in geothermal wells (Vasini et al., 2017). This version contains an analytical solution of heat transfer between the wellbore and surrounding rock formations following the approach by Ramey (1962) which includes both the steady-state heat conduction through the components of well completion (casings and cement sheaths) and the radial transient heat conduction within the rock formation. Ramey’s analytical function for heat exchange between wellbore and formation is substituted by the Chiu and Thakur (1991) time function.
phase which is computed neglecting the pressure effects. This approach is acceptable for low partial pressures of CO2, at which its enthalpy in the gas phase can be evaluated with an ideal gas approach. When the partial pressure becomes high, the CO2 gas enthalpy depends also on pressure, and the pressure effect is simply transmitted to the enthalpy of dissolved CO2 if the heat of solution is computed neglecting the pressure. Koschel et al. (2006) showed that the enthalpy of CO2 dissolved in the aqueous phase is almost constant with respect to pressure, indicating that the effect of pressure on the heat of solution is almost compensating for that on gaseous CO2 enthalpy. Thus, in addition of the previous approach a new one can be optionally chosen to calculate the enthalpy of dissolved CO2 as a function of temperature only. - Under sub-critical conditions, phase transitions taking place during TOUGH2 numerical simulations are occasionally affected by numerical instabilities which, among other issues, may be related to the difficulties to interpolate the fluid properties at the interface between grid blocks. The latter problem has been studied by O’Sullivan et al. (2013, 2014) for the EOS modules commonly employed for geothermal reservoir simulation, including EOS2. In order to obtain a more stable interpolation of phase density across a phase transition, they implemented the calculation of the density of non-existing phases. The same approach is now optionally available in EOS2H. In addition, the same option performs also the calculation of enthalpy and viscosity of the non-existing phase, in a way similar to that implemented into iTOUGH2 (Finsterle, 2015).
3. T2well-EOS2H verification 3.1. Verification of reservoir flow against iTOUGH2-EOS1sc The initial verification of EOS2H has been conducted by using the TOUGH2 RFP sample problem (Pruess et al., 1999) involving the injection and production at constant rate of 30 kg/s using a 5-spot well pattern. Taking advantage of the symmetry, the grid is limited to 1/8 of the well pattern and is built parallel to the diagonal connecting injection and production wells. EOS2H has successfully reproduced (Saipem, 2016) (a) the original RFP sample problem results simulated for pure water with EOS1 and EOS2 with reservoir T of 300 °C, (b) the EOS2 results at 300 °C and initial PCO2 of 0.5 and 5 MPa, and (c) the iTOUGH2-EOS1sc results for pure water with initial reservoir
In addition, EOS2H can optionally simulate a traced water component in complete analogy to EOS1 (Pruess et al., 1999). The option can be employed to trace reinjected water, such as the separated brine and/ 3
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 3. P vs T conditions simulated halfway between injection and production wells: comparison between results obtained with iTOUGH2-EOS1sc (symbols and lines) and TOUGH2-EOS2H (lines). Modified after Magnusdottir and Finsterle (2015). IAPWS-IF97 values computed by EOS1sc are overlapped by those computed by EOS2H.
correlation to the other has not shown any numerical problems in the simulated cases: with the maximum time step of 36.5 days, no time step reductions were encountered, nonlinearities were solved with a reasonable number of Newtonian iterations.
conditions of 1100 °C and 90 MPa, 1500 °C and 150 MPa, and 1200 °C and 50 MPa (Magnusdottir and Finsterle, 2015). Comparison of EOS2H results and those simulated with iTOUGH2-EOS1sc are shown in Fig. 3 for the grid element located halfway between the injector and producer. The somewhat unrealistic initial reservoir supercritical conditions (with temperatures above the melting point of commonly encountered rocks) were chosen to replicate the previous iTOUGH2-EOS1sc results, and to show that extrapolated correlations for Regions 2 and 5 do work properly. iTOUGH2-EOS1sc results are correctly replicated by EOS2H. The PT path of selected grid elements of the 3 simulations is reported in Fig. 4 on the PT thermodynamic diagram of pure water, on which the limits of IAPWS-IF97 regions are also drawn. The PT path of production and reinjection wells is shown, as well as that of the element halfway between the two wells. In the T = 1500 °C case, the production well block remains in the extrapolated field of Region 5, while the injection well block goes from the extrapolated field of Region 2 inside Region 3 and then back into the extrapolated field of Region 2. In the T = 1200 °C case, the intermediate element and the injection well block go from the extrapolated field of Region 5 inside Region 2, while the production well block remains in the same Region. In the T = 1100 °C case, the intermediate element and the production well block remain in Region 5, while the injection one enters the extrapolated field of Region 5 to terminate inside Region 2. The simulated cases have involved Regions 2 and 5 and their extrapolated fields above 100 and 50 MPa, respectively, which are the upper limits of the IAPWS-IF97 correlations. The switch from one
3.2. Verification of reservoir flow against AUTOUGH and STAR-HOTH2O Additional simulations have been performed to compare EOS2H with the AUTOUGH results of Croucher and O’Sullivan (2008) who reproduced previous simulations performed with the STAR-HOTH2O simulator (Pritchett, 1994) by Yano and Ishido (1998). A 1D radial model is used with outer radius of 12,039 m, thickness of 100 m, discretized with 40 elements with a radius increment ratio of 1.3. The model features and rock properties listed by Croucher and O’Sullivan (2008) have been used. Two series of simulations were performed with a constant production rate of 15.7 kg/s for 300 h at: A.) reservoir P of 30 MPa and T from 200 to 500 °C; and B.) reservoir T of 450 °C and P from 25 to 35 MPa. The comparison of EOS2H results and those of STAR (black lines) and AUTOUGH (black symbols) are shown in Figs. 5 and 6. Case A at 400 °C (Fig. 5) cannot be modeled with EOS2H as reservoir conditions corresponding to Region 3 are encountered. For the other cases a good reproduction of AUTOUGH results is shown, while the departure from STAR results at early times is due to a difference in near well 4
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 4. P&T thermodynamic diagram of pure water: evolution of reservoir conditions at the injection and production well locations and at the halfway point for the 3 simulated cases.
analysis to the simulated drawdown results. Simulations have been performed starting with 23 MPa initial reservoir pressure and a CO2 partial pressure of 2 MPa, and with initial temperatures of 350, 370, 400, 450 and 500 °C. The partial pressure of water is then equal to 21 MPa, which is slightly lower than the water saturation pressure at 370 °C equal to 21.0434 MPa. Then, the simulation at 370 °C has been performed by assigning two-phase initial conditions with a gas phase saturation of 0.5; the simulation at 350 °C refers to single-liquid conditions; those at 400, 450 and 500 °C refer to supercritical steam-like conditions. In all cases the well is operated at a constant mass extraction rate of 15.7 kg/s for 300 h (1.08 × 106 s). The resulting pressures as a function of time in the well block
permeability introduced by Croucher and O’Sullivan (2008). Fig. 6 also shows a good reproduction of AUTOUGH and STAR results, except for early time P computed with STAR by Yano and Ishido (1998). The comparison confirms that the IAPWS-IF97 correlations have been correctly included into TOUGH2, and that no adverse numerical problems are encountered using EOS2H. 3.3. Code demonstration against a simple drawdown solution The reservoir model used in Section 3.2 above has been employed with different initial reservoir conditions to simulate short production tests of 300 h in order to apply a customary semi-log pressure transient
Fig. 5. Wellbore P vs time: reservoir P 30 MPa. EOS2H vs AUTOUGH (symbols) and STAR (lines) results (modified after Croucher and O’Sullivan, 2008). 5
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 6. P change vs time: reservoir T 450 °C. EOS2H vs AUTOUGH (symbols) and STAR (lines) results (modified after Croucher and O’Sullivan, 2008).
Fig. 7. P = 23 MPa and CO2 partial pressure of 2 MPa. Well block pressure is shown as function of production time on a linear scale.
evaluation of reservoir hydraulic transmissivity (m3) with the following equation:
element are shown in Figs. 7 and 8, on linear and logarithmic scales, respectively. The pressure drawdown at the well block increases with increasing temperature, except for the 370 °C case, as the two-phase conditions lead to a greater compressibility than the single-liquid case at 350 °C. The fluid density decreases with increasing temperatures, and then the drawdown is greater when the extraction rate is held constant. The semi-log plot in Fig. 8 highlights linear trends at late time corresponding to the radial flow dominated period which can be used for the evaluation of permeability-thickness product kh in analogy to the pressure data recorded downhole during a drawdown pressure transient. In Fig. 8 the data at late time regressed with a straight line are shown as function of the logarithm of time. For each case the regression straight line has been determined with the evaluation of the slope m (Pa/cycle) which can be used for the
kh = 2.3026
Wμ 4πρ m
(1)
where: k permeability, m2 h layer thickness, m W mass rate, kg/s μ dynamic viscosity, Pa s ρ density, kg/m3 Tables 1 and 2 list the slope values of the regression straight line, the permeability k evaluated with Eq. 1 using the fluid viscosity and density at i) initial reservoir conditions and ii) in the well block at the 6
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 8. P = 23 MPa and CO2 partial pressure of 2 MPa. Well block pressure is shown as function of logarithm of production time.
changes more at increasing reservoir temperature, and a temperature decline on the order of 10 °C is also simulated as shown in Fig. 9. As viscosity and density of steam and CO2 mixture depend on both pressure and temperature, they change more due to reservoir exploitation. Eq. 1 is strictly valid only for a low compressibility fluid with constant thermo-physical properties. Thus, the difference between the computed permeability value and that used for the numerical simulations is due to the non-negligible variations of supercritical mixture properties during reservoir pressure drawdown.
Table 1 Slope of regression straight line and kh product evaluated with viscosity and density values corresponding to initial reservoir conditions for single-phase cases.
m (Pa/cycle) μ (Pa s) q (m3/s) ρ (kg/m3) kh (m3) k (m2) Error k (%)
350°C
400°C
450°C
500°C
3.42E+05 7.15E-05 2.55E-02 6.17E+02 9.76E-13 9.76E-15 −2.37
7.35E+05 2.67E-05 1.25E-01 1.26E+02 8.32E-13 8.32E-15 −16.81
1.05E+06 2.84E-05 1.59E-01 9.88E+01 7.87E-13 7.87E-15 −21.27
1.40E+06 3.04E-05 1.84E-01 8.55E+01 7.29E-13 7.29E-15 −27.10
3.4. Verification of steady state wellbore flow against T2Well-EWASG and PROFILI simulators T2Well-EOS2H has been compared to T2Well-EWASG by simulating the steady-state wellbore flow in two geothermal wells discharging fluids with high contents of CO2; details are given by Vasini et al. (2018). The simulations are related to:
Table 2 Slope of regression straight line and kh product evaluated with viscosity and density values corresponding to final well block reservoir conditions for singlephase cases.
m (Pa/cycle) μ (Pa s) q (m3/s) ρ (kg/m3) kh /m3) k (m2) Error k (%)
350°C
400°C
450°C
500°C
3.42E+05 7.03E-05 2.58E-02 6.08E+02 9.73E-13 9.73E-15 −2.68
7.35E+05 2.54E-05 1.60E-01 9.81E+01 1.02E-12 1.02E-14 +1.58
1.05E+06 2.73E-05 2.18E-01 7.22E+01 1.03E-12 1.03E-14 +3.40
1.40E+06 2.93E-05 2.72E-01 5.78E+01 1.04E-12 1.04E-14 +3.79
a) Well KD-13, Kizildere field, Turkey, discharging a low NaCl brine (1000 ppm) with a high content of dissolved CO2 (20,000 ppm); data published by James (1975); b) Well W2 from an unspecified field managed by ENEL, discharging a brine with remarkable NaCl (9600 ppm) and CO2 (30,000 ppm) contents; data published by Barelli et al. (1982). The results of the comparison are shown in Fig. 10, where field measurements are presented together with T2Well-EWASG and T2WellEOS2H results. Even though slight differences are present, which are most likely related to a different wellbore gridding and to the need to adjust the CO2 concentration due to the absence of NaCl effects, the results show that T2Well-EOS2H is able to reproduce both the experimental steady-state wellbore flow data and the corresponding T2WellEWASG simulation results of Vasini et al. (2018). T2Well-EOS2H is also compared to a steady-state wellbore flow simulation performed with PROFILI (Saipem, 2017) based on conditions of Well IDDP-1, Krafla field, Iceland (Einarsson et al., 2015). The numerical wellbore flow simulator PROFILI (Battistelli, 2010) simulates the steady-state flow under production or injection of H2ONaCl-CO2 fluid mixtures under single-phase liquid, two-phase and single-phase gas conditions (including supercritical gas-like conditions)
end of 300 h of discharge, respectively. While the permeability for the 350 °C single-phase liquid is evaluated with a comparable error for the different fluid parameters (−2.37 and −2.68 %, respectively), for the cases with supercritical temperature, the error on permeability is much higher, between −16 and −27 % when using the initial reservoir conditions, with respect to that obtained using the final well block fluid parameters, which is between 1.58 and 3.79 %. For the single-phase liquid case at 350 °C, no phase transitions are experienced and then both fluid viscosity and density remain almost constant as the pressure transient occurs at almost isothermal conditions and liquid fluid properties do not change considerably with pressure. For supercritical steam-like conditions, the well block pressure 7
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 9. Well block P&T path for the 5 simulated cases: reservoir pressure of 23 MPa and CO2 partial pressure of 2 MPa.
4. Modeling of coupled wellbore-reservoir flow at well sector scale
in vertical or directional wellbores with variable diameters and hole surface roughness. The numerical solution of mass, energy and momentum balance equations is based on the work by Barelli et al. (1982). For two-phase conditions, pressure changes due to friction, acceleration and gravity can be evaluated using a homogeneous approach based on the DIF-3 correlations by CISE (Lombardi and Ceresa, 1978; Bonfanti et al., 1979) or the CESNEF-2 correlations (Lombardi and Carsana, 1992). Well IDDP-1 is completed with a 9-5/8” casing down to 1950 m and with a 9-5/8” slotted liner down to 2070 m. The wellbore is here assumed to be vertical even though the production tests were performed using the final wellbore path resulting from two side-tracks. The PROFILI simulation is performed from well-head to bottom-hole. The rate is 7 kg/s with wellhead P (WHP) of 14.2 MPa, wellhead T (WHT) of 440 °C, and 1000 ppm CO2. Heat transfer is computed following Ramey (1962) with an assumed average global heat exchange coefficient U of 20 W/(m2 °C). The T2Well-EOS2H simulation is performed by fixing flowing bottom conditions equal to those computed by PROFILI and extracting 7 kg/s from the well-head until almost steady-state conditions are reached after 30 days. Fig. 11 shows good agreement between flowing P&T computed by PROFILI and T2Well-EOS2H.
A preliminary coupled wellbore-reservoir flow simulation was previously presented by Battistelli et al. (2018) based on the results of production tests performed at well IDDP-1, Krafla, Iceland (here referred as 2018 results). The simulations did not intend to accurately replicate the IDDP-1 production test results as testing of the well was complex and performed in 5 stages of quite different duration from March 2010 to June 2012. The production testing assemblage was changed according to encountered operating conditions with some long shut-in periods between production stages. In addition, Einarsson et al. (2015) provide the WHP, WHT and wellhead enthalpy (WHE) history as functions of time but not the related discharge rate history. Thus, because of the complexity of IDDP-1 testing operations and the lack of some important data, the simulations did not intend to accurately reproduce the observed well behavior, but just consider IDDP-1 published results as a realistic framework to be used for the modeling of coupled wellbore-reservoir flow. First, the well output curve (OC) reconstructed by Einarsson et al. (2015) was reproduced using PROFILI by adjusting the static reservoir P&T and using a parabolic P draw-down as a function of mass rate with
Fig. 10. Simulated P&T vs depth: (a) well KD13; (b) well W2. T2Well-EOS2H vs T2Well-EWASG. 8
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 11. Simulated P&T vs depth for well IDDP-1: T2Well-EOS2H vs PROFILI comparison.
while both WHT and WHE are slightly higher due to the higher overall heat exchange coefficient of well completion with respect to the average value used previously. The OC is then reproduced with T2Well-EOS2H by simulating in full transient conditions well discharge at increasing rates from 5 to 50 kg/s with constant production steps lasting 6 days each as in the new PROFILI steady state simulation. The reservoir is described using the same 1D radial model of 1 km radius and 100 m thickness (from 2000 to 2100 m depth) previously employed. The results in terms of WHP and WHT for 2018 and 2019-A simulations are shown in Fig. 13; they are obtained by trial and error using permeabilities of 8.45 × 10−14 m2 and 8.70 × 10−14 m2, respectively. Porosity has been set at 0.10 and no skin effects have been included. WHP and WHT are plotted vs discharge time, while the steady-state values computed using PROFILI are shown as circles at the end of each production step. WHP was previously (Battistelli et al., 2018) slightly underestimated by T2Well-EOS2H for rates of 15–30 kg/s and overestimated at the highest rate of 50 kg/s, when the fluid velocity simulated by PROFILI reached 209 m/s at the well-head, being only 51 m/s at bottom-hole. At simulated wellhead conditions the fluid velocity was still much lower than the speed of sound of pure water which is about 646 m/s. The friction gradient increased from 2.50 × 103 Pa/m at bottom-hole to 7.05 × 103 Pa/m at the wellhead. At 45 kg/s the fluid velocity increased from bottom to well-head from 40.6 m/s to only 67.0 m/s. Thus, the increment from 45 to 50 kg/s produced an important increment of fluid velocity which was nevertheless still well below chocked flow conditions. The new 2019-A simulation shows similar trends for WHP and WHT, with WHP lower than PROFILI at rates lower than 40 kg/s, and higher at rates above 40 kg/s. The increment of production step duration from 3 to 6 days induced a non-negligible inflow of supercritical fluid from the lateral outer boundary at the highest rates. In order to have an almost infinite acting reservoir, a new grid was built with the lateral boundary set at 5 km, and with the same discretization as the previous grid in the first 1000 m. The BHP and BHT computed with T2Well-EOS2H (called 2019B) are also shown in Fig. 13. They were obtained by increasing the reservoir permeability to 9.30 × 10−14 m2. Apart from that, the results of the simulation 2019-A are very close to those of simulation 2019-B. The latter results are used in the subsequent comparison with the PROFILI 2019 results.
viscous and turbulent coefficients of 0.5 bar/(kg/s) and 0.0175 bar/ (kg/s)2, respectively. The global heat exchange coefficient of well completion, U, was assumed constant at 10 W/(m2 °C) to compute the heat exchange with surrounding rock formations following Ramey (1962). The heat conductivity and diffusivity of the surrounding rock formations were assumed constant with depth and set to 2.25 W/(m °C) and 9.3E-7 m2/s, respectively. The IDDP-1 production casing is a Hydrill 563, 9-5/8″ 53.5 lb/ft pipe, while the slotted liner is a 9-5/8″ 47 lb/ft pipe with Buttress connections (Pálsson et al., 2014). The roughness of casing and liner was assumed to be 1.0E-5 m and 1.0E4 m, respectively. The output curve was computed at increasing rates from 5 to 50 kg/s, assuming a constant average production time of 30 days for each rate. While the bottom hole P (BHP) changes with rate, the bottom hole T (BHT) was held constant for each production step, and a PCO2 of 7.2 kPa corresponding to 1000 ppm was assumed. Initial static P&T values were calibrated at 15.5 MPa and 510 °C. The OC was then reproduced with T2Well-EOS2H by simulating well discharge at increasing rates from 5 to 50 kg/s with constant production steps lasting only 3 days each, instead of 30 days. The reservoir was described using a 1D radial model of 1 km radius and 100 m thickness (from 2000 to 2100 m depth) with homogeneous, isotropic rock properties, initial P&T and fluid composition equal to that calibrated using the PROFILI simulation. The radial grid includes the wellbore block plus 100 radial elements with the lateral boundary at 1 km kept at its initial conditions. All the parameters related to the heat exchange were the same as those used for PROFILI simulations. For the new PROFILI simulation described below (referred to as 2019 results), the global heat exchange coefficient of well completion, U, is computed following Willhite (1967) and using the final well design given by Pálsson et al. (2014) at 10.25 W/(m2 °C) down to 254 m, 12.19 W/(m2 °C) down to 785 m, and 26.86 W/(m2 °C) down to 1900 m. The output curve is computed at increasing rates from 5 to 50 kg/s, with 10 steps of 6 days each and the production time for each step is varied accordingly. The same parameters describing the downhole deliverability curve already calibrated in Battistelli et al. (2018) were used. Fig. 12 shows the 2018 and 2019 results of PROFILI simulations with the OC and WHT on the left, and the WHE on the right. Changes made in the calculation of heat exchange had minor effects on the OC which is almost unchanged with respect to the previous simulation, 9
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 12. IDDP-1 output curve: mass rate and WHT vs WHP (a) and mass rate and WHE vs WHP (b) simulated using PROFILI (modified after Einarsson et al., 2015). “Calculated mass Flow” and “Characteristic Curve” are the experimental rates and interpreted output curve, respectively, of Einarsson et al. (2015).
high rates. In fact, water specific heat at 15 MPa and T from 410 to 510 °C changes from 3.91 to 2.85 kJ/(kg °C), while it is substantially higher for subcritical liquid water as it ranges between 4.15 and 5.47 kJ/(kg °C) from 50 to 300 °C. Geothermal wells discharging supercritical fluids are likely to experience high heat losses in the relatively cold upper well section which turns into a relevant T reduction of discharged fluid due to its low heat capacity. Among the reasons for the WHT differences between PROFILI and T2Well we can first mention that while PROFILI uses a steady-state approach, T2Well computes a transient solution even when the heat transfer is computed with the analytical method by Ramey (1962). Moreover, the PROFILI simulations have been performed by assuming a constant BHT, as the simulator just computes the wellbore flow. On the other hand, T2Well-EOS2H computes the BHT as a result of reservoir flow by solving both mass and heat balance equations. Fig. 16 shows
A direct comparison between the 2019 output curve simulated by PROFILI and T2Well-EOS2H (2019-B) is shown in Fig. 14, where both mass rate and WHT are plotted versus WHP. Fig. 15 shows WHP and WHT versus mass rate. The T2Well-EOS2H WHP is lower and higher than PROFILI above and below 40 kg/s, respectively. The maximum discrepancy is obtained at the highest rate of 50 kg/s where the WHP computed by T2Well-EOS2H is significantly higher. Fig. 15 shows how the WHT increases considerably with rate up to about 15 kg/s mostly because the higher rates imply smaller heat losses per unit mass. Then the WHT almost levels out up to about 40 kg/s. It finally declines at the highest rates. Differences between the two sets of results is within ± 5 °C, apart for the 45 kg/s results. Temperature changes for wells discharging supercritical water seem to be relatively sensitive to enthalpy variations, as those related to the heat exchange with surrounding formation and kinetic energy losses at
Fig. 13. IDDP-1 output curve: WHP and WHT vs time simulated using T2Well-EOS2H. 2018 (a), 2019-A (b) and 2019-B (c) simulation results compared with corresponding steady-state PROFILI results plotted at the end of each step. 10
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 14. IDDP-1 output curve: Mass rate and WHT vs WHP simulated using T2Well-EOS2H, compared with values computed by PROFILI and the experimental ‘Characteristic Curve’ by Einarsson et al. (2015).
Fig. 15. IDDP-1 output curve: WHP and WHT vs mass rate simulated using T2Well-EOS2H, compared with values computed by PROFILI.
could be a better assumption for the PROFILI wellbore flow steady-state simulations. Then, the OC was simulated again using PROFILI but assuming a constant BHE equal to the initial reservoir value. The computed BHT as function of rate is shown in Fig. 18 which suggests that with the assumption of a constant BHE the BHT are lower than those actually computed by T2Well-EOS2H. Fig. 19 shows the effect of the constant BHE assumption on the PROFILI output curve: the computed WHT becomes much lower than the previous 2019 results, while the effects on the WHP are significant only at the highest rate of 50 kg/s. In practice, because of the strong dependency of supercritical steam enthalpy on both P&T, the common assumption of assigning the BH conditions in wellbore flow simulation
the BHP and BHT as function of rate for both PROFILI and T2Well. The BHT computed by T2Well declines at increasing rates, becoming 57.7 °C lower at the highest rate of 50 kg/s. As can be expected, the enthalpy at bottom hole has also a different trend for the two simulations. Starting with the same value of 3330 kJ/ kg at initial reservoir conditions, PROFILI computes a final BHE of 3414 kJ/kg greater than the 3359 kJ/kg of T2Well, because of the PROFILI assumption of constant BHT. Fig. 17 shows that the WHE increases at increasing rates up to the maximum simulated rate of 50 kg/s. The increment is higher up to 10−15 kg/s. Then the WHE increment is reduced as the rate increases. The BHE history of T2Well suggests checking if a constant BHE
Fig. 16. IDDP-1 downhole deliverability curve: BHP and BHT vs mass rate simulated using T2Well-EOS2H, compared with values used by PROFILI. 11
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Fig. 17. IDDP-1 output curve: WHP and WHE vs mass rate simulated using T2Well-EOS2H, compared with values computed by PROFILI.
Fig. 18. IDDP-1 downhole deliverability curve: BHP and BHT vs mass rate simulated using T2Well-EOS2H, (OC 2019-B) compared with BHT values computed by PROFILI when assuming constant BHE.
EOS2H is capable of accurately modeling coupled wellbore-reservoir flow processes under supercritical steam-like conditions.
for wells producing from liquid dominated reservoirs do not simply apply when modeling the flow of supercritical fluids. The coupled wellbore-reservoir flow modeling seems to be able to supply more reliable results than a wellbore flow simulator. While the simulated temperatures might be considerably affected by the accuracy of heat balance calculations and related assumptions, computed flowing pressures seem to be less affected. These preliminary simulations of coupled wellbore-reservoir flow at a well sector scale show promising results, suggesting that T2Well-
5. Conclusions The existing EOS2 module of the TOUGH2 V.2.0 reservoir simulator for H2O-CO2 mixtures has been improved by including the capability to model the thermodynamic conditions of supercritical steam-like reservoirs using the IAPWS-IF97 correlations within Regions 2 and 5. The
Fig. 19. IDDP-1 output curve: Mass rate and WHT vs WHP simulated using T2Well-EOS2H and the experimental output curve (Einarsson et al., 2015), compared with values computed by PROFILI using both constant BHT and BHE assumptions. 12
Geothermics 86 (2020) 101793
A. Battistelli, et al.
Validation, Writing - original draft. Stefan Finsterle: Methodology, Writing - review & editing. Marica Marcolini: Writing - review & editing, Project administration. Lehua Pan: Methodology, Writing review & editing.
calculation of CO2 properties at high T and P has been improved by including the use of PR EOS for density and enthalpy departure as well as the LBC method for the dynamic viscosity. The improved module, called EOS2H, has been linked to T2Well for the modeling of coupled wellbore-reservoir flow under fully transient conditions in geothermal reservoirs under both sub-critical and supercritical steam-like conditions (Regions 2 and 5 of the IAPWS-IF97 correlations). T2Well-EOS2H accuracy and numerical performances have been successfully verified by comparison with:
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments
- TOUGH2-EOS2 V.2.0 reservoir simulations at sub-critical conditions; - iTOUGH2-EOS1sc reservoir simulations at supercritical conditions without CO2; - STAR-HOTH2O and AUTOUGH reservoir simulations at supercritical conditions; - T2Well-EWASG and PROFILI wellbore flow simulations.
This work was supported by Saipem SpA under the R&D project "Modeling of flow in geothermal wells", Job P5000Q-06. The management of Saipem SpA is acknowledged for the permission to publish this paper. References
The verified code has then been used to model the output curve of a well producing from a supercritical steam-like geothermal reservoir. Wellbore and reservoir characteristics are based on those of well IDDP1 drilled in the Krafla field, Iceland. The output curve of well IDDP-1 inferred from discharge tests has been replicated using the PROFILI steady-state numerical wellbore simulator with capabilities to simulate steam-like supercritical flow conditions, by calibrating the coefficients of the downhole deliverability curve. Then the output curve has been reproduced with T2Well-EOS2H using a simple model with a 1D vertical grid for the wellbore and a 1D radial grid for the supercritical reservoir. The numerical model calibration has been limited to a simple trial and error approach by modifying reservoir permeability. A more rigorous approach was deemed to be not necessary because several of the wellbore and reservoir features used were not available and then were just assumed based on educated guesses. T2Well-EOS2H simulated output curve shows differences with PROFILI results which are likely to be due to:
Barelli, A., Corsi, R., Del Pizzo, G., Scali, C., 1982. A two-phase flow model for geothermal wells in the presence of non-condensable gas. Geothermics 11 (3), 175–191. Battistelli, A., Calore, C., Pruess, K., 1997. The simulator TOUGH2/EWASG for modelling geothermal reservoirs with brines and a non-condensible gas. Geothermics 26 (4), 437–464. Battistelli, A., 2010. PROFILI Wellbore Flow Simulator: From Geothermal Wells to GHG and Acid Gas Injection Wells. GeoThermExpo 2010, Ferrara September 21-22, 2010. Battistelli, A., Swenson, D., Alcott, A., 2017. Improved PetraSim-TOUGH2 Capabilities for the Simulation of Geothermal Reservoirs. 42nd Stanford Geoth. Work, Feb. 13-15, 2017, Stanford, CA. Battistelli, A., Finsterle, S., Pan, L., 2018. Modeling of coupled wellbore-reservoir flow for supercritical geothermal systems in steam-like conditions. In: TOUGH Symposium 2018. Berkeley, CA. Bonfanti, F., Ceresa, I., Lombardi, C., 1979. Two-phase pressure drops in the low flowrate region. Energia Nucleare 26 (10), 481–492. Brikowski, T.H., 2001. Modeling supercritical systems with TOUGH2: preliminary results using the EOS1SC equation of state module. In: Proc. 26th Work. Geoth. Res. Eng., Stanford University. Stanford, CA. Chase Jr, M.W., 1998. NIST-JANAF themochemical tables. J. Phys. Chem. Ref. Data, Monograph 9, Fourth Edition. pp. 1–1951. Chiu, K., Thakur, S.C., 1991. Modeling of wellbore heat losses in directional wells under changing injection conditions. In: SPE 22870. Presented at SPE Annual Technical Conference and Exhibition. October, 1991. Dallas, Texas. pp. 517–528. Cox, B.L., Pruess, K., 1990. Numerical experiments on convective heat transfer in watersaturated porous media at near-critical conditions. Transp. Porous Media 5 (3), 299–323. Croucher, A.E., O’Sullivan, M.J., 2008. Application of the computer code TOUGH2 to the simulation of supercritical conditions in geothermal systems. Geothermics 37, 622–634. Duan, Z.H., Sun, R., 2003. An improved model calculating CO2 solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar. Chem. Geol. 193, 257–271. Einarsson, K.K., Sveinsson, E., Ingason, K., Kristjánsson, V., Hólmgeirsson, S., 2015. Discharge testing of magma Well IDDP-1. In: Proc. WGC 2015. Melbourne, Australia.. Ellis, A.J., Golding, R.M., 1963. The solubility of carbon dioxide above 100°C in water and in sodium chloride solutions. Am. J. Sci. 261, 47–60. Finsterle, S., Commer, M., Edmiston, J., Jung, Y., Kowalsky, M.B., Pau, G.S.H., Wainwright, H., Zhang, Y., 2017. iTOUGH2: a multiphysics simulation-optimization framework for analyzing subsurface systems. Comput. Geosci. 108, 8–20. Finsterle, S., 2015. Enhancements to the TOUGH2 Simulator Implemented in iTOUGH2, Report LBNL-7016E. Lawrence Berkeley National Laboratory, Berkeley, CA. IAPWS, 2007. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam. Switzerland. IAPWS, 2008. Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance. Berlin, Germany. James, R., 1975. Gas content of a hot-water reservoir estimated from downhole pressure and temperature measurements. In: Berkeley, CA. Second United Nations Symposium Vol. 3. pp. 1689–1691. Kissling, W.M., 1995. Extending MULKOM to super-critical temperatures and pressures. In: Proc. WGC 1995. Firenze, Italy. pp. 1687–1690. Kissling, W.M., 2005. Transport of three-phase hyper-saline brines in porous media: theory and code implementation. Transp. Porous Media 61, 25–44. Koschel, D., Coxam, J.-Y., Rodier, L., Majer, V., 2006. Enthalpy and solubility data of CO2 in water and NaCl(aq) at conditions of interest for geological sequestration. Fluid Phase Equilib. 247, 107–120. Lohrentz, J., Bray, B.G., Clark, C.R., 1964. Calculating viscosities of reservoir fluids from their compositions. SPE paper 915. J. Pet. Technol. Altern. Fuels 1171–1176. Lombardi, C., Ceresa, I., 1978. A generalized pressure drop correlation in two-phase flow. Energia Nucleare 25 (4), 181–198. Lombardi, C., Carsana, C.G., 1992. A dimensionless pressure drop correlation for two-
i) the different numerical solutions implemented, with steady-state wellbore flow for PROFILI, and full transient wellbore-reservoir flow for T2Well-EOS2H, and ii) the use of a downhole deliverability curve for PROFILI with the assumption of fixing bottom hole flowing conditions as either constant BHT or constant BHE, while variable BHT and BHE are actually computed by T2Well-EOS2H. The simulations of wellbore flow in the presence of supercritical fluids suggest that an accurate modeling of heat transfer toward the surrounding rock formation might be more important than for conventional sub-critical geothermal wells. In fact, geothermal wells discharging supercritical fluids are likely to experience large heat losses in the relatively cold upper well section which turns into a significant temperature reduction of the discharged fluid due to its low heat capacity. This aspect requires additional efforts in order to improve the present modeling capabilities, even though using a full numerical approach for the modeling of radial heat conduction would probably be already able to improve the reliability of the results. The preliminary simulations of coupled wellbore-reservoir flow presented above show promising results, suggesting that T2Well-EOS2H is capable of accurately modeling the coupled wellbore-reservoir flow under supercritical steam-like conditions as those encountered near magmatic heat sources and within the deep roots of already exploited geothermal fields. CRediT authorship contribution statement Alfredo Battistelli: Conceptualization, Methodology, Software, 13
Geothermics 86 (2020) 101793
A. Battistelli, et al. phase mixtures flowing up in vertical duct covering wide parameter ranger. Heat Pipe Sci. Technol. Int. J. 10 (1-2), 125–141. Magnusdottir, L., Finsterle, S., 2015. An iTOUGH2 equation-of-state module for modeling supercritical conditions in geothermal reservoirs. Geothermics 57, 8–17. McKibbin, R., McNabb, A., 1999. Deep hydrothermal systems: mathematical modeling of hot dense brines containing noncondensible gases. J. Porous Media 2 (1), 107–126. Montegrossi, G., Pasqua, C., Battistelli, A., Mwawongo, G., Ofwona, C., 2015. 3D natural State model of the menengai geothermal system, Kenya. In: Proc. WGC 2015. Melbourne, Australia, 19–25 April 2015. Morin, A., Lövfall, B.T., Meese, E., 2016. Simulation of Supercritical Water Flow in the Venelle 2 Well. (Accessed 24 June 2019). http://www.descramble-h2020.eu/ images/download/publications/GeoThArt.pdf. NIST, 2007. Standard Reference Database 23. NIST Thermophysical Properties of Hydrocarbon Mixtures – Program REFPROP – Version 8.0.Standard Reference Database 23. NIST Thermophysical Properties of Hydrocarbon Mixtures – Program REFPROP – Version 8.0. O’Sullivan, J., Croucher, A., Yeh, A., O’Sullivan, M.J., 2013. Improved convergence for air-water and CO2-water TOUGH2 simulations. In: Proceedings 35th New Zealand Geoth. Work. Rotorua, New Zealand. O’Sullivan, J., Croucher, A., Yeh, A., O’Sullivan, M.J., 2014. Further improvements in the convergence of TOUGH2 simulations. In: Oñate, E., Oliver, J., Huerta, A. (Eds.), 11th World Congress on Computational Mechanics (WCCM XI). O’Sullivan, J., Kipyego, E., Croucher, A., Ofwona, C., O’Sullivan, M.J., 2015. A supercritical model of the menengai geothermal system. In: Proc. WGC 2015. Melbourne, Australia, 19–25 April 2015. O’Sullivan, J., O’Sullivan, M.J., Croucher, A., 2016. Improvements to the AUTOUGH2 supercritical simulator with extension to the air-water equation-of-State. GRC Transactions 40, 921–930. O’Sullivan, M.J., Bodvarsson, G.S., Pruess, K., Blakeley, M.R., 1985. Fluid and heat flow in
gas-rich geothermal reservoirs. SPE J. 25, 215–226. Pálsson, B., Hólmgeirsson, S., Guðmundsson, Á., Bóasson, H.Á., Ingason, K., Sverrison, H., Thórhallsson, S., 2014. Drilling of the well IDDP-1. Geothermics 49, 23–30. Pan, L., Oldenburg, C.M., 2014. T2Well - An integrated wellbore–reservoir simulator. Comp. Geosci. 65, 46–55. Peng, D.-Y., Robinson, D.B., 1976. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 15, 59–64. Pritchett, J.W., 1994. HOTH2O A Description of H2O Properties to one Kilobar and 800°C for Use With the STAR Geothermal Reservoir Simulator. SSS-TR-94-14730, S-Cubed, La Jolla. Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User’s Guide, Version 2.0, Report LBNL-43134. Lawrence Berkeley National Laboratory, Berkeley, CA. Ramey Jr, H.J., 1962. Wellbore heat transmission. In: 36Th Ann. Fall Meet. of SPE, Oct. 811, 1961. Dallas. Saipem, 2016. Numerical Modeling of Hydrothermal Systems with Supercritical Steam Conditions. Module EOS2H, Report SPC_94060, Job P5000Q-06_03 (unpublished). . Saipem, 2017. Simulation of Production Tests With TOUGH2-T2Well, Report SPC_94070, P5000Q-06_03 (unpublished). . Spycher, N., Pruess, K., 2010. A phase-partitioning model for CO2–brine mixtures at elevated temperatures and pressures: application to CO2-Enhanced geothermal systems. Transp. Porous Media 82, 173–196. Vasini, E.M., Battistelli, A., Berry, P., Bonduà, S., Bortolotti, V., Cormio, C., Pan, L., 2018. Interpretation of production tests in geothermal wells with T2Well-EWASG. Geothermics 73, 158–167. Willhite, G.P., 1967. Over-all heat transfer coefficients in steam and hot water injection wells. J. Pet. Technol. Altern. Fuels 607–615 Paper SPE 1449. Yano, Y., Ishido, T., 1998. Numerical investigation of production behavior of deep geothermal reservoirs at supercritical conditions. Geothermics 27 (5/6), 705–721.
14