Calculation of temperature profile in injection wells

Injection wells have long been an essential asset in enhanced oil recovery, wastewater disposal and carbon dioxide sequestration in petroleum industries. The temperature profile of fluid flow in the injection well is one of the main parameters of interest for petroleum engineers to determine optimum injection conditions and wellbore completion design especially in thermal injection projects and deep wells. In this study, the calculation involved in determining the temperature profile along the depth of wellbore has been revised to be newer and more robust via solving governing wellbore equations. The wellbore is segmented into discrete counterparts for it to be solved simultaneously in terms of mass, momentum and energy balance via wellbore governing equations. Five injection cases from the literatures, incompressible and compressible fluid flows, were used to confirm that the procedure is reproducible in terms of its behaviour, which is similar to field data. The new results acquired from the new procedure are in good agreement with field data collected with a maximum absolute error less than 3 °C.


List of symbols a
Geothermal gradient (°C/m) A Area (m 2 ) A c Coefficient matrix A tubi Internal area of the tubing (m 2 ) A r Coefficient (m) A z Internal area of the tubing at level z (m 2 ) A z+Δz Internal area of the tubing at level z + Δz (m 2 ) b Surface geothermal temperature (°C) C an Specific heat of the fluid in the annulus (J/(kg °C)) C e Heat capacity of the earth (J/(kg °C)) C j Joule Thomson coefficient (°C/Pa) C p Specific heat of the fluid at constant pressure (J/ (kg °C)) e Internal energy per unit mass (J/kg) E t Energy content within the element at time t (W) E t+Δt Energy content within the element at time t + Δt (W) e z Internal energy per unit mass at level z (J/kg) e z+Δz Internal energy per unit mass at level z + Δz (J/kg) F Momentum losses due to friction (N) f Moody friction factor g Acceleration of gravity (m/s 2 ) G r Grashof number h Enthalpy per unit mass (J/kg) h c Heat transfer coefficient for convection (W/ (m 2 °C)) h f Film coefficient (W/(m 2 °C)) h r Heat transfer coefficient for radiation (W/(m 2 °C)) h z Enthalpy per unit mass at level z (J/kg) h z+Δz Enthalpy per unit mass at level z + Δz (J/kg) i Number of cells in the radius direction k Kinetic energy per unit mass (J/kg) k a Thermal conductivity of the fluid in the annulus (W/(m °C)) k c Equivalent thermal conductivity of the fluid in the annulus (W/(m °C)) k cas Thermal conductivity of the casing (W/(m °C)) k cem Thermal conductivity of the cement (W/(m °C)) k e Thermal conductivity of the earth (W/(m °C)) k ins Thermal conductivity of the insulation (W/(m °C)) k tub Thermal conductivity of the tubing (W/(m °C))

Introduction
Injection wells have long been an essential asset in enhanced oil recovery (EOR), wastewater disposal and carbon dioxide sequestration in petroleum industries (Hasan et al. 2002;Moradi 2013;Hamdi et al. 2014). The temperature of fluid flow in wellbore is one of the main parameters of interest for petroleum engineers. In-depth understanding of the pressure and temperature profiles along the depth of the well is a requirement for the appropriate design of well. The performance of hydrocarbon reservoirs can only be gauged with precise determination of downhole pressure and temperature. In addition to that to prevent damaging the formation by injecting above threshold pressure, extensive knowledge of the bottomhole pressure is useful. Although temperature can be measured by bottomhole gauge, the possibility of downhole gauge failure increases over a long period of time. Thus, the ability to calculate downhole parameters from surface injection parameters would be of great convenience (Paterson et al. 2008).
Currently, there are several accurate software packages that are available to calculate the temperature profile of flowing fluid along the depth of wellbore based on computational fluid dynamics (CFD) solutions (Fluent 2011). However, CFD solutions are rather slow and need a high capacity machine to run it; there is also the difficulty in model building for inexperienced users (Apak 2006). Also there are other types of software packages (Wellflo 2001; VFPi 2011) available for the determination of pressure and temperature profiles in the wellbore based on Ramey's model (Ramey Jr 1962). The running speed of these packages is fast but Ramey's model has been developed based on some assumptions that are not suitable for fluid flow at near critical point or deep injection wells (Messer et al. 1974;Alves et al. 1992;Yasunami et al. 2010). In addition to that there are some users who have reported difficulties in defining fluid thermodynamic properties in some of these commercial wellbore simulators.
In order to overcome these issues, this study is objectively conducted to develop a rapid and reliable procedure to determine pressure and temperature profiles free from the aforementioned limitations.

Methodology
Simultaneous solving of equations governing the wellbore in terms of momentum, mass and energy balance has to be carried out for the calculation of temperature profile (Herrera et al. 1978;Fontanilla and Aziz 1982;Hasan et al. 2002;Paterson et al. 2008). The main components involved in the systems include tubing, insulation, casing and formation. Other elements involved are the fluid flowing through the tubing, annulus between the insulation and casing as well as the cement surrounding the casing (Moradi 2013). To begin finding resolution for the equations governing the wellbore, the system shown in Fig. 1 is discretized into N segments in the vertical direction to consider various fluid thermodynamic properties, overall coefficient of heat transfer, rate of heat transfer and heterogeneity of layers around the wellbore along the wellbore depth. The optimum value for length of segments is calculated by performing a sensitivity analysis (Yasunami et al. 2010). It is assumed that all variables within a segment remain constant. As shown in Fig. 2, the pressure and temperature of every segment are (Livescu et al. 2010): (1) The density (ρ f ), viscosity (µ f ) and velocity (u f ) of the fluid in the segment can be calculated by knowing T f , P f and using a fluid thermodynamic properties table (Peng and Robinson 1976; National Institute of Standards and Technology 2011). For the calculation of fluid properties in the thermodynamic module, data published by the National Institute of Standards and Technology (NIST) is used for pure materials, and Peng-Robinson equation of state (Peng and Robinson 1976) is applied for calculating mixture properties.
Comparatively, the heat transferred along the wellbore is faster than the heat transferred in the layers surrounding the wellbore and the formation attributed to the small wellbore radius. Besides, the small wellbore radius also contributed to the rate of heat transfer reaching a steady state much sooner than in the formation as heat is transferred under unsteady state (Ramey Jr 1962;Fontanilla and Aziz 1982). Therefore, steady-state rate of heat transfer is made to solve the wellbore governing equations, while unsteady-state rate of heat transfer is made for the formation governing equation, without introducing significant errors.

Mass balance equation
The calculation of velocities at level z and z + ∆z is done using mass balance equation. In a steady-state condition and a given volume (Fig. 2), the flow of fluid within the tubing in accordance with the conservation of mass law is (Kreith et al. 2010): where ṁ is defined by: The combination of Eqs. (3) and (4) calculates the velocities at level z and z + ∆z as follows: The radius of the tubing is constant (Ali 1981;Fontanilla and Aziz 1982). The Internal area of the tubing equals:

Momentum balance equation
In a steady-state condition, momentum equation is the basic governing equation utilized for the calculation of its pressure gradient. The total forces acting on the fluid element is equivalent to the momentum change in the fluid; in the wellbore, the forces present are pressure, friction and gravity (Hasan et al. 2002;Pan et al. 2007;Paterson et al. 2008;Livescu et al. 2010). For a given control volume (Fig. 2): The term F is the loss in momentum as a result of friction, and it is depicted by (Hasan et al. 2002;Pan et al. 2007;Paterson et al. 2008;Livescu et al. 2010): By substituting (7) and (9) in (8), simplifying, rearranging and dividing both sides of (9) by A tubi , there will be: For calculating pressure by (10), the Moody friction factor (f) should be defined is usually an expression of Reynolds number: Generally, empirical correlations in computer calculations are used to determine Moody friction factor. Laminar flow is defined when Reynolds number (Re) is < 2400, friction factor is inversely related to Reynolds number [40].
Turbulent flow is defined when Reynolds number (Re) ≥ 2400. In this condition, the friction factor is dependent on both Reynolds number and pipe roughness (Chen 1979;Pan et al. 2007). In a turbulent flow condition, the empirical correlation presented by Chen (1979) for determining f is: where

Energy Balance Equations
The difference in temperatures between wellbore fluid flow and strata around the wellbore leads to an exchange of energy. The energy balance equation for every segment must be solved to calculate heat transfer rate in each element of the system as shown in Fig. 2 (Ramey Jr 1962;Hasan et al. 2002Hasan et al. , 2010. In a steady-state system, the general energy balance is (Ramey Jr 1962;Hasan et al. 2002;Kreith et al. 2010):  (15) for the given control volume (Fig. 2), there is: The definition of enthalpy gives (Kreith et al. 2010): Rearrangement of (17) and replacement of Pv + u by h yield: In (18), heat generation rate is only accounted for by a single term on the right due to the loss in flow friction and the rate of conductive heat transfer between wellbore and surrounding formation. Nevertheless, the rate of heat generated by flow friction loss is so minute that it is negligible (Pan et al. 2007;Paterson et al. 2008).
By rearranging (18), the following equation is obtained: In a steady-state condition, the rate of heat flowing through a wellbore of the segment is written as (Hasan et al. 2010

Heat flow rate through the earth and wellbore/formation interface
The surrounding of the wellbore system that acts as a heat source or sink is known as the formation. The transfer of heat between wellbore fluid and the surrounding earth is attributed to their temperature difference. For modelling heat flow in the earth, first, the temperature profile has to be modelled in the formation. Heat transfer in the formation around the wellbore is dictated by the conduction mechanism instead of the convection mechanism (Ramey Jr 1962;Hasan and Kabir 1991;Moradi 2013;Moradi and Awang 2013). The Laplace transform was used by Hasan and Kabir (1991) to model the temperature profile of the formation. This was done following an approach taken by Van Everdingen and Hurst (1949) for pressure transients using a comparable set of equations. Thus, the definition of dimensionless temperature, T D , is: where They assumed surrounding layers around the wellbore are homogeneous with constant heat flow rate at the wellbore/formation interface. It is discovered that the following algebraic expressions for dimensionless temperature, T D , in terms of dimensionless time, t D , have quite accurately represented the solution.

Fluid temperature
The assumption that temperature, pressure and flow velocity are constant in a cross section of inside the tubing (T f , P f and v f ) is made due to the small ratio between tubing diameter and length. Thus, the flow is considered being one-dimensional (Hasan et al. 2002;Yasunami et al. 2010 (Yasunami et al. 2010).

Tubing temperature
High conductivity metal is used to make the tubing [e.g. k of steel is 43.275 W/(m °C)]; hence, the temperature distribution is considered negligible and T tubi = T tubo (Ali 1981).

Annulus temperature
There is a temperature profile across the annulus. Normally, an average temperature of the outer surface of the insulation and inner surface of the casing is assumed as annulus temperature, T an (Pourafshary et al. 2009;Hasan et al. 2010):

Casing temperature
High conductivity metal is used to make the casing [e.g. k of steel is 43.275 W/(m °C)]; hence, the temperature distribution is considered negligible and T casi = T caso (Yasunami et al. 2010); so by neglecting the thermal resistance of the casing:

Calculation steps of the new procedure
As stated in the objectives, a procedure is presented to predict the fluid flow pressure and temperature profiles in the injection well by using equations derived in the previous sections. Input data needed to run the new procedure include P, T and ṁ f of injection flow at wellhead, r tubi , r tubo , r inso , r casi , r caso , r cemo , k ins , k cem , k e , e , C e , T s , tubi , C an , an , k a , inso , and casi . The well is discretized into N segments and a procedure as shown below is use in all segments: Step 1 Assign a random value for T at level z + ∆z. The temperature at level z can be used for initial guess.
Step 2 Assign a random value for P at level z + ∆z. The pressure at level z can be used for initial guess.
Step 3 Using a thermodynamic module determines the density, enthalpy and viscosity of fluid at level z and z + ∆z.
Step 4 Determine the average properties of the segment at P f and T f by (1) and (2), respectively. Then, calculate ρ f , µ f and v f at T f and P f in the thermodynamic module.
Step 5 Resolve momentum equation and calculate P at level z + ∆z by (10).
Step 6 Compare the calculated P at level z + ∆z from step 5 with the assumed P from step 2. If (calculated P − assumed P) < 0.001 (kPa), the procedure has converged and vice versa. When non-convergent of the procedure occurs, return to step 2, assumed P = calculated P and repeat the procedure until same values are obtained.
Step 7 Assign a random value for the ∆Q/∆z (for initial guess ∆Q/∆z = 0).
Step 9 Equate T casi to the geothermal temperature, T e .
Step 10 Calculate the value T ins by (31).
Step 14 Compare the calculate T casi from step 13 with step 9, if ABS(old T casi − new T casi ) < 0.01 (°C), the procedure is assumed to have converged and go to step 17 else go to step 15.
Step 15 Determine the corresponding ∆Q/∆z from (20) based on the calculated T cemo from step 12 and go to Step 10 and repeat until convergent of the procedure.
(34) T cemo = r tubo * U to * T D * T f + k e * T e r tubo * U to * T D + k e Step 16 Determine the specific enthalpy at level z + ∆z by (19).
Step 17 Determine the T z+∆z at P z+∆z and h z+∆z in developed thermodynamic module.
Step 18 Compare the calculated T at level z + ∆z from step 17 with the assumed T from step 1, if ABS(calculated T − assumed T) < 0.01 (°C), convergent of the procedure is assumed. If there is a non-convergent of procedure, return to step 1, assumed T = calculated T and repeat until convergent of procedure is obtained. Figure 3 represents the flowchart of the calculation procedure.

Validation
Validation is a necessity in ensuring the validity of the mathematical formulation, solution techniques, and program coding. The validity of the new procedure is examined against three water injection cases, one carbon dioxide injection case and one steam injection case to confirm that the procedure is reproducible in terms of its behaviour, and in agreement with collected field data. A computer code in Visual C#.net environment has been programmed in order to easily implement the new procedure shown in Fig. 3 for various cases.
Case 1 Figure 4 compares results of the new procedure against the field data taken from Nowak (1953). Water at a surface temperature of 28.33 °C was injected at the rate of 0.0016 m 3 /s through 0.1778 m casing diameter for three years. There, the calculated temperature exceeded the measured temperatures, but the difference is small. The maximum discrepancy over the depth is 2.39 °C.
Case 2 For further verification on the performance of the procedure presented, newly calculated results using the procedure are compared with those presented by Squier et al. (1962) Case 3 As a third case, the computer code was examined with cold water injection, data from Ramey Jr (1962). Water was injected at 0.0088 m 3 /s through 0.162 m casing diameter for a period of 75 days. Injection temperature was 14.7222 °C at the wellhead. The results from the new procedure agree quite well with Ramey's results as shown in Fig. 6. Comparison of the results predicted by the new procedure with Ramey's data showed that the maximum discrepancy over the depth is 0.81 °C.
Case 4 For the confirmation of the validity of the proposed procedure for compressible fluid, a comparison was made with the data presented in the study of Cronshaw and Bolling (1982). At the time of the survey, carbon dioxide at a surface temperature of 13.39 °C was injected at a rate of 163.871 m 3 /s. Figure 7 revealed a fairly good match between bottomhole temperatures of carbon dioxide after 7 days of injection as reported by Cronshaw and Bolling with that predicted by the new procedure. Cronshaw and Bolling reported the bottomhole temperature as 25.22 °C. The present procedure predicted a value of 26.36 °C.
Case 5 For further validation of the model, results calculated using the new procedure are compared with those presented by Satter (1965). Superheated steam was injected into a vertical well at a mass flow rate of 0.63 kg/s at 537.78 °C and 3.45 MPa. Figure 8 showed a good agreement between the calculated results and Satter's results after 3.65 days injection. Comparison of the results predicted by the new procedure with the results reported by Sattar showed an average absolute per cent relative deviation of 0.15.

Application of the new procedure
Input data needed to run the new procedure are listed as the following: • Surface parameters: P, T and ṁ f of injection flow at wellhead. I. Radius: r tubi , r tubo , r inso , r casi , r caso and r cemo . II. Thermal properties: k ins , k cem , tubi , C an , an , k a , inso casi .
• Formation properties: T s , k e , ρ e and C e .
The main outputs of the new procedure are pressure and temperatures profiles along the depth of wellbore.
Beside those, the new procedure calculates the below items: • Geothermal temperature versus depth. • Heat transfer through each facility of the well completion. • Heat transfer through layers around the wellbore. • Kinetic energy, potential energy. • Contribution of changing the kinetic energy and potential energy for building the temperature profile. • Contribution of changing the hydrostatic gradient, acceleration gradient, and frictional gradient for building the pressure profile.
The new procedure that was proposed in this study is useful for the following applications: • Comparing measured data with calculated data This could be for one of several purposes, such as evaluating "match parameters" which are difficult or impossible to measure, pipe roughness, or determining if a well is behaving the way it is expected to (i.e. to detect faulty components). • Monitoring work such as predicting bottomhole pressure from measured surface pressure and flow rate.

Fig. 4
Comparison of temperature distribution with well depth from Nowak (1953) and predicted by the new procedure • Conducting a design work where it is required to calculate the pressure and temperature drop in the wellbore such as: to determine the best diameter of the tubing and injection flow rate.
The coupled response from fluid movement in the wellbore to the reservoir is neglected that is the main limitation of the new workflow. This is important during early stage of production or injection as the fluid flows under unsteady conditions specially during well testing. To handle this issue, the new producer must be coupled to a reservoir simulator. Also it needs to add unsteady-state formulations to calculate the energy balance.

Conclusions
The main output of this study is a non-isothermal wellbore simulator. This new simulator discretizes the wellbore to several segments and considers various thermodynamic properties and overall heat transfer coefficient for every segment. The ability of the new simulator to calculate its own overall heat transfer coefficient holds substantial benefit over commercial software packages. The validity of the new procedure and computer code has been examined in various scenarios against the results from the literatures and they agreed quite well with each other.   Cronshaw and Boling (1982) and predicted by the new procedure Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.