Thermal Profiles in Water InjectionWells: Reduction of the Systematic Error in the Flow Measurement during the Transient Regime

Preprint

Article

Thermal Profiles in Water InjectionWells: Reduction of the Systematic Error in the Flow Measurement during the Transient Regime

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

Submitted:

16 September 2023

Posted:

20 September 2023

You are already at the latest version

Alerts

Abstract

Every simplification has the potential to add a systematic error to the measured value. Sometimes, the addition does not render the measurement useless. In the case of oil well flow measurements based on the thermal properties of the fluid in the well, the addition of a systematic error of up to 20% can occur if the variation in the internal energy of the fluid and, most importantly, the transient nature of heat transfer in the completion are not taken into account. This error is a mathematical consequence of the simplifications chosen. To work within the first hours of good operation, given the need for immediate application in a prototype water injection well developed and installed at the UFRN, this article presents an analytical solution for calculating the flow rate by solving the differential equations obtained by applying the principles of mass and heat transfer. As a result, at the cost of greater complexity, the systematic error drops to values of less than 1% in the first two hours of operation of the well, as seen throughout this document.

Keywords:

Subject: Engineering - Electrical and Electronic Engineering

1. Introduction

In the search for increasing or maintaining oil production in mature fields, several techniques for reenergizing reservoirs and readjusting the physical properties of oil (viscosity, surface tension, ...) have been developed throughout the industry’s history. Among them, the water injection technique was one of the first to be studied and disseminated, being the most used today [1].

Wells intended for the injection of fluids, called injection wells, are distinguished from wells intended for oil extraction, called producing wells. The fluid used is called injection fluid. They generally have the same type of completion (tubing, annulus, casing, and cementing), as represented in the schematic in Figure 1.

Measuring the injected flow rate is crucial for the excellent performance of any injection method. With this information, it is possible to detect deviations from the designed flow rates caused by failures in the mechanical flow regulators and the damage existing to the formation [2]. When adopting multiple injection zones, in which fluid is simultaneously injected into several sections along the injection well, the quality of this measurement is more critical, given the risk of greater deviation from the designed flow.

Flow measurement using thermal profiles depends on the existence of a temperature gradient and, consequently, two thermal profiles. In the case of injection wells, the thermal profile established in the injection fluid and the geothermal profile, established naturally in the formation, have been used. The first has a defined location, in the middle of the pipe. The second has an undefined location, tending to infinity, in accordance with the progression of the thermal disturbance. With the knowledge of these two profiles, the flow rate will be completely defined, if the physical properties (conductivity, thermal capacity, etc.), the operating time, the initial conditions, and the well boundary conditions are known.

Given the impossibility of physically accessing all points of the well and reservoir to collect the physical properties and initial condition, the following conceptual simplifications have been adopted in the literature since the first studies [3]:

(a): formation with infinite extension;
(b): linear, time-invariant geothermal profile;
(c): unidirectional and constant injection flow rate;
(d): constant fluid temperature at the surface;
(e): zero vertical heat flow;
(f): homogeneous, isotropic, and time-invariant physical properties;
(g): initial temperatures equal to the geothermal profile.

As a result of the simplifications listed, which have been adopted in whole or in large part in all of the technical texts referenced in this article, the solutions developed have low accuracy in the first moments of an injection well’s operation. In particular, the quasi-static methods, which include the original Ramey [4] solution, require weeks to obtain a meaningful reading. The analytical methods studied, Wu & Pruess [3] and Assmann [5], because they do not take into account the transient nature of heat transfer in the completion, also suffer from low accuracy in flow measurement in the first few moments.

2. Preliminaries

Consider an idle injection well in thermal equilibrium with the formation, as shown in Figure 1. Assume that a new injection process starts at

t = 0

. Also, assume that the injection fluid has no phase change, is incompressible, and is colder than the reservoir. As the fluid moves through the pipe, it gains thermal energy, and its temperature increases, creating a thermal profile along the pipe. The established profile is directly related to the flow rate in the pipe.

To calculate the flow rate for

t > 0

, after the well operation’s start, consider the principle’s energy and mass conservation application in the Control Volume highlighted in Figure 2. Assume that vertical heat flow can be neglected and that all thermal properties of the well, formation, and fluid, as well as all initial and boundary conditions, are known. Also, assume that the flow rate is constant and the flow pattern is turbulent. Due to this last characteristic, the fluid can be addressed as an agglomerated system, in which the transience of heat transfer in the fluid can be disregarded, given the slowness of transfer mechanisms in your neighborhood. Consider also that the thermal properties are isotropic, homogeneous, and invariant over time and that there is no heat generation due to the fluid viscosity.

As a result, there is [3,6]:

\begin{matrix} \dot{m} ≅ \frac{2 π r_{1} U (z, t)}{c} [T_{g} (z) - T_{f} (z, t)] / \frac{\partial T_{f} (z, t)}{\partial z} - π {(r_{1})}^{2} ρ \frac{\partial T_{f} (z, t)}{\partial t} / \frac{\partial T_{f} (z, t)}{\partial z} \end{matrix}

(1)

where

U (z, t)

is the Global Heat Transfer Coefficient, with

r = r_{1}

as a reference. It encompasses all thermal resistance to heat propagation, from convection on the pipe’s internal surface, at

r = r_{1}

, through pipe conduction and natural convection in the annulus, to the radial limit of thermal disturbance caused by transporting the fluid at a colder temperature.

The calculation of the flow rate,

\dot{m}

, as can be deduced from Equation (1), depends on knowing the value of

U (z, t)

, which in turn, as it results from the thermal disturbance carried out, depends on understanding the convection coefficient on the internal surface of the pipe, which depends on the flow rate you want to measure. In other words, calculating the flow rate will depend on solving the following system of non-linear equations:

P r = \frac{μ c}{k}

(2)

R e = \frac{ρ υ_{m e d} (2 r_{1})}{μ}

(3)

f = {[0.79 l n (R e) - 1.64]}^{- 2}

(4)

h = \frac{(f / 8) (R e - 1000) P r}{1 + 12.7 {(f / 8)}^{0.5} ({P r}^{2 / 3} - 1)} \frac{k}{2 r_{1}}

(5)

\begin{matrix} \tilde{U} ≅ U (\tilde{z}, t) \end{matrix}

(6)

\dot{m} ≅ \frac{2 π r_{1} \tilde{U}}{c} \frac{\int_{z_{1}}^{z_{2}} [T_{g} (z) - T_{f} (z, t)] d z}{T_{f} (z_{2}, t) - T_{f} (z_{1}, t)} - π {(r_{1})}^{2} ρ \frac{\int_{z_{1}}^{z_{2}} [\partial T_{f} (z, t) / \partial t] d z}{T_{f} (z_{2}, t) - T_{f} (z_{1}, t)}

(7)

where:

P r

is the Prandtl number,

R e

is the Reynolds number, f is the friction factor, h is the convection coefficient,

U (\tilde{z}, t)

is the global average heat transfer coefficient in the interval [

z_{1}

z_{2}

] and

\dot{m}

is the desired mass flow rate. The variables

μ

, c, k, and

ρ

are, respectively, the absolute viscosity, the specific heat capacity, the thermal conductivity, and the specific mass, all properties relative to the fluid in flow.

The equation (4) corresponds to Petukhov’s first explicit equation for calculating the friction factor in smooth tubes. It is valid for

3 * 10^{3} < R e < 5 * 10^{6}

([7]).

The equation (5) corresponds to the Gnielinski equation for calculating the convection coefficient in forced turbulent flow inside pipes. It is valid for

0.5 \leq P r \leq 2 * 10^{3}

and

3 * 10^{3} < R e < 5 * 10^{6}

[7].

The system resulting from equations (2) to (7) can be solved, among other methods in the literature, by executing the algorithm shown in Figure 3.

The approach chosen for calculating

\tilde{U}

will define its applicability and the intensity of the corresponding systematic error. This article will present two approaches: a) the quasi-static approach, in which it is assumed that the thermal processes are already slow enough to be considered in thermal equilibrium, and b) the purely analytical approach, in which the solution is sought by solving the differential heat transfer equations.

3. Quasi-Static Approach

Given the small radius of the completion compared to the formation and the long operating times of the injector and producer wells, the quasi-static solutions disregard all thermal transience in the fluid and the completion (tubing, annulus, casing, and cementing). As a result, still taking Figure 2 as a reference, the Global Heat Transfer Coefficient and the Mass Flow will be given by:

\begin{matrix} \tilde{U} ≅ {[\frac{1}{h} + \frac{r_{1}}{k_{t u b}} l n (\frac{r_{2}}{r_{1}}) + \frac{r_{1}}{r_{2}} \frac{1}{h_{c o m b}} + \frac{r_{1}}{k_{r e v}} l n (\frac{r_{4}}{r_{3}}) + \frac{r_{1}}{k_{c i m}} l n (\frac{r_{5}}{r_{4}}) + \frac{r_{1} f (t)}{k_{f o r}}]}^{- 1} \end{matrix}

(8)

\dot{m} ≅ \frac{2 π r_{1} \tilde{U}}{c} \frac{\int_{z_{1}}^{z_{2}} [T_{g} (z) - T_{f} (z, t)] d z}{T_{f} (z_{2}, t) - T_{f} (z_{1}, t)}

(9)

where

f (t)

is the Transient Function for heat propagation in the formation up to the thermal perturbation limit [4]. This function encompasses the entire time dependence of heat propagation in the formation. It is the primary variable to be determined in this methodology. It is defined by:

\begin{matrix} f (t) = \frac{2 π k_{f o r} [T (r_{5}, z, t) - T_{g} (z)] Δ z}{Δ \dot{Q} (r_{5}, z, t)} \end{matrix}

(10)

In Equation (10),

T (r_{5}, z, t)

is the temperature at the interface with the formation,

T_{g} (z)

is the geothermal temperature at depth z and

Δ \dot{Q} (r_{5}, z, t)

is the heat rate transferred to the formation by the depth element

Δ z

[4].

In 1962, Ramey proposed the following approximation for times longer than a week [4]:

f (t) ≅ - l n \frac{r_{5}}{2 \sqrt{α t}} - 0.290

(11)

with

α

being the thermal diffusivity of the formation and t being the time measured in days.

In 1967, Mattews & Russell, assuming that the fluid inside the pipe could be considered a linear geometric heat source (line source), given the small radius of the pipe compared to the radius of the formation, came up with the following approximation for the transient function [9]:

f (τ) ≅ \frac{1}{2} \int_{\frac{1}{4 τ}}^{\infty} \frac{e^{- u}}{u} d u

(12)

with

τ

corresponding to the dimensionless Fourier time for formation, i.e:

τ = \frac{k_{f o r} t}{ρ_{f o r} c_{f o r} {r_{5}}^{2}}

(13)

In 1994, Hasan & Kabir proposed the following approximation for the transient function referring to the condition of constant heat flow at the interface with the formation, in

r_{1}

[8]:

τ \leq 1.5 \to f (τ) = 1.1281 \sqrt{τ} (1 - 0.3 \sqrt{τ})

(14)

τ > 1.5 \to f (τ) = [0.4063 + 0.5 l n (τ)] (1 + \frac{0.6}{τ})

(15)

In 2004, Hagoort proposed an approximate function for the constant temperature boundary condition at the internal interface of the formation. This is [10]:

τ < 1 \to f (τ) = \frac{1}{\frac{1}{π τ} + 0.5 - 0.2 \sqrt{τ / π}}

(16)

1 < τ < 10^{5} \to f (τ) = - 0.0012 {[l n (τ)]}^{3} + 0.0249 {[l n (τ)]}^{2} + 0.3083 l n (τ) + 1.0504

(17)

τ > 10^{5} \to f (τ) = 0.5 l n (\frac{τ}{0.447})

(18)

In general, the transient functions for boundary conditions of the type “constant temperature at

r_{1}

”, “constant flow at

r_{1}

” and “fluid with constant temperature and radiation boundary condition at

r_{1}

” are given, respectively, by:

\begin{matrix} f (τ) = 1 / L_{s \to τ}^{- 1} {\frac{K_{1} (\sqrt{s})}{\sqrt{s} K_{0} (\sqrt{s})}} \end{matrix}

(19)

\begin{matrix} f (τ) = L_{s \to τ}^{- 1} {\frac{K_{0} (\sqrt{s})}{s \sqrt{s} K_{1} (\sqrt{s})}} \end{matrix}

(20)

\begin{matrix} f (τ) = L_{s \to τ}^{- 1} {\frac{1}{s} \frac{K_{0} (\sqrt{s})}{K_{0} (\sqrt{s}) + \frac{\sqrt{s}}{β_{f o r}} K_{1} (\sqrt{s})}} / L_{s \to τ}^{- 1} {\frac{1}{\sqrt{s}} \frac{K_{1} (\sqrt{s})}{K_{0} (\sqrt{s}) + \frac{\sqrt{s}}{β_{f o r}} K_{1} (\sqrt{s})}} \end{matrix}

(21)

where

K_{0}

and

K_{1}

are the Bessel functions of the second kind and order 0 and 1, respectively, and

β_{f o r}

is the Biot number for the formation [7], i.e:

β_{f o r} = \frac{r_{1} U_{f o r}}{k_{f o r}}

(22)

with

U_{f o r}

being the Global Heat Transfer Coefficient up to the interface with the formation, given by:

\begin{matrix} U_{f o r} = {[\frac{1}{h} + \frac{r_{1}}{k_{t u b}} l n (\frac{r_{2}}{r_{1}}) + \frac{r_{1}}{r_{2}} \frac{1}{h_{c o m b}} + \frac{r_{1}}{k_{r e v}} l n (\frac{r_{4}}{r_{3}}) + \frac{r_{1}}{k_{c i m}} l n (\frac{r_{5}}{r_{4}})]}^{- 1} \end{matrix}

(23)

The derivation of the listed equations, i.e., the solution of the system of differential equations for the mentioned boundary conditions, can be found in the work of [6].

4. Analitical Approach

Wu & Pruess in 1990 [3], Assmann [5] in 1993, and Haggort in 2004 [10] presented analytical solutions for calculating the temperature of the fluid inside the pipe. Using the concept of the global heat transfer coefficient to synthesize the thermal resistances up to the formation, they all disregarded the transience in completion. In addition, the work by Wu & Pruess presented data on the relationship between vertical and radial gradients, justifying the assumption of zero vertical heat flow.

Although the results obtained by Wu & Pruess, Assmann, and Hagoort, among others, show good results after the well has been in operation for some time, disregarding the transience of heat transfer in the completion, as well as the variation in the internal energy of the fluid (second part of Equation (7), makes it impossible to obtain better results for times around the transit time of the fluid in the pipe, which corresponds to the time required for the fluid to travel all the way through the pipe to the depth under analysis.

Consider the simplified completion shown in Figure 4, where the annulus has been ignored to illustrate the methodology. The layers of the well have been numbered to facilitate the naming of properties and variables.

As we have seen so far, calculating the flow rate depends on knowing the value of the overall heat transfer coefficient, which in turn depends on the convection coefficient, the fluid temperature, the geothermal temperature, and the temperature of the internal pipe interface. Consider the following dimensionless temperatures for the fluid and the internal pipe interface:

θ_{f} (z, t) = \frac{T_{f} (z, t) - T_{g} (z)}{T_{g} (z)}

(24)

θ_{1} (z, t) = \frac{T (r_{1}, z, t) - T_{g} (z)}{T_{g} (z)}

(25)

Since the heat flow by convection at the internal interface of the pipe must be equal to the heat flow at the same interface due to the use of the global heat coefficient, we have:

U (z, t) = h \frac{T_{f} (z, t) - T (r_{1}, z, t)}{T_{f} (z, t) - T_{g} (z)}

(26)

and, consequently:

U (z, t) = h - \frac{θ_{1} (z, t)}{θ_{f} (z, t)} h

(27)

Considering that the vertical thermal gradient, and consequently the heat flow in the same direction, is only significant when compared to the intensity of the radial flow at the thermal disturbance front, which occurs during the fluid’s transit time in the pipe, it will be disregarded in this analysis [3]. Thus, excluding the vertical heat flow, assuming that the flow regime inside the pipe is turbulent and applying the principles of conservation of energy and mass, we arrive at the system of partial differential equations described in Table 1 [6].

The deduction of the equations and the solution of the resulting system can be found in Appendix A of the work by Lima [6]. The estimates for the temperature of the fluid and the temperature of the internal surface of the pipe are given by:

\begin{matrix} θ_{f} (w, τ) = \frac{1}{T_{g} (w)} [- a_{w} L_{s \to τ}^{- 1} {\frac{1 - e^{- g (s) w}}{s g (s)}} + b_{w} θ_{f 0} L_{s \to τ}^{- 1} {\frac{e^{- g (s) w}}{s}}] \end{matrix}

(28)

and

\begin{matrix} θ_{1} (w, τ) = \frac{1}{T_{g} (w)} [- a_{w} L_{s \to τ}^{- 1} {{\bar{\bar{θ}}}_{1} (w, s) \frac{1 - e^{- g (s) w}}{s g (s)}} + b_{w} θ_{f 0} L_{s \to τ}^{- 1} {{\bar{\bar{θ}}}_{1} (w, s) \frac{e^{- g (s) w}}{s}}] \end{matrix}

(29)

where:

\begin{matrix} g (s) = 2 α_{1} β_{1} [1 - {\bar{\bar{θ}}}_{1} (w, s)] + s \end{matrix}

(30)

\begin{matrix} {\bar{\bar{θ}}}_{1} (w, s) = {\bar{c}}_{11} I_{o} (\sqrt{s}) - {\bar{c}}_{12} K_{o} (\sqrt{s}) \end{matrix}

(31)

such that:

\begin{matrix} {\bar{θ}}_{1} (w, s) = {\bar{\bar{θ}}}_{1} (w, s) {\bar{θ}}_{f} (w, s) \end{matrix}

(32)

The function

g (s)

is used to simplify the notation. The functions

I_{0} (\sqrt{s})

and

K_{0} (\sqrt{s})

are known in the literature as modified Bessel functions of the first and second kind, respectively, both of zero order. The variables

{\bar{c}}_{11}

and

{\bar{c}}_{12}

are functions in the Laplace domain, independent of the radial coordinate, r, and related to the mean 1 of the completion under study, taken as a reference in the process of dimensioning the variables. Their values are given by the system of linear functions in the Laplace domain, which is a consequence of the boundary conditions of the problem. See [6] for more information.

In equations (28) and (29),

τ

is the dimensionless time, w is the dimensionless vertical coordinate,

θ_{f 0}

is the dimensionless temperature of the fluid at

w = 0

a_{w}

is the geothermal gradient related to the dimensionless coordinate w and

b_{w}

is the linear coefficient of the dimensionless geothermal temperature. The values of

τ

, w,

a_{w}

and

b_{w}

are given, respectively, by:

\begin{matrix} τ = \frac{k_{1}}{ρ_{1} c_{1} {r_{1}}^{2}} t \end{matrix}

(33)

\begin{matrix} w = \frac{k_{1}}{v_{m e d} ρ_{1} c_{1} {r_{1}}^{2}} z \end{matrix}

(34)

\begin{matrix} a_{w} = \frac{v_{m e d} ρ_{1} c_{1} {r_{1}}^{2}}{k_{1}} a_{z} \end{matrix}

(35)

\begin{matrix} b_{w} = T_{g} (0) \end{matrix}

(36)

The dimensionless variables represented by the Equations (24), (25), (33) and (34) are consequences of the algebraic manipulations that attempt to reduce the mathematical complexity of the system of partial differential equations. Among them, the dimensionless time variable (

τ

) is known in the literature as dimensionless Fourier time.

Thus, in addition to equation (27), the global heat coefficient can be rewritten as:

U (z, t) = h - \frac{- a_{w} L_{s \to τ}^{- 1} {{\bar{\bar{θ}}}_{1} (w, s) \frac{1 - e^{- g (s) w}}{s g (s)}} + b_{w} θ_{f 0} L_{s \to τ}^{- 1} {{\bar{\bar{θ}}}_{1} (w, s) \frac{e^{- g (s) w}}{s}}}{- a_{w} L_{s \to τ}^{- 1} {\frac{1 - e^{- g (s) w}}{s g (s)}} + b_{w} θ_{f 0} L_{s \to τ}^{- 1} {\frac{e^{- g (s) w}}{s}}} h

(37)

The inverse Laplace transforms of Equation (37) can be calculated numerically using the Gaver-Stehfest algorithm [11], which was adopted in this article, or other methods available in the literature.

As seen, it is possible to use a standard script for both approaches by continuing to use the global heat transfer coefficient as the basis for calculating the flow rate. However, unlike the previous quasi-static approach, the calculation of the global coefficient given by Equation (37) also takes into account the thermal transience in the layers of the well.

5. Computer Simulation

Given the very specific operating conditions, with data collection and storage in the first moments of injection, in a reservoir with original initial conditions, the possibility of the existence of real, unmotivated data, from the operation of a well in production under the conditions listed has proved to be nil up to the date of completion of this article, making it considerably difficult to evaluate the proposed analytical solution.

As a countermeasure to the lack of data for validation, a computerized well simulator was designed and developed to generate the temporal evolution of the thermal profile of the injection fluid along the pipe.

Due to its widespread use in the scientific community, the finite difference method was used as a numerical tool in the simulator code development. In this method, the central idea is to replace the differential equations with algebraic equations, exchanging the derivatives for difference approximations and applying the resulting equations to each subdivision of the problem’s physical space, referred to in this article as the Control Volume (CV). The result is an algebraic equation for each subdivision, making up a system of linear equations [7].

The code was developed in Phyton, using the Jupyter Notebook editor, version

5.4 . 0

, from Anaconda Navigator, version

1.7 . 0 .

See the code at [6].

6. Results

Two computer simulations were carried out to evaluate the proposed method. The first simulation aimed to verify the alignment of the proposed solution with the classical quasi-static solutions proposed by Ramey and his successors. To this end, a schematic of a well without completion was adopted, aligning it with the reference model used to calculate the transient functions. The second simulation aimed to verify the proposed analytical solution in a well with simplified completion (with tubing and cementing) and compare its performance with quasi-static solutions.

6.1. Simulation 1: Well without Completion

In this first analysis, the data generated in the simulation of a well without completion is processed, as shown in Figure 5. This configuration is a good representation of the prototype well installed at UFRN. Although the prototype has a tubing, its high thermal conductivity means that its presence does not make a noticeable difference in the fluid temperature evolution, as will be seen below.

The well shown in Figure 5 cannot be used as a real-world model. However, by processing the fluid temperature evolution for this configuration, it is possible to evaluate the contribution of the transient functions to the flow measurement, since these functions were derived from a similar schematic of the well. Furthermore, considering that water injection wells operate for months, sometimes years, heat conduction in the formation is the dominant thermal process in flow measurement. Under these conditions, after the first two weeks of operation, any thermal transient in the completion that could influence the flow inference becomes insignificant [4]. In this way, a simulation without completion provides indications and trends for the behavior of the flow to be measured.

6.2. Simulation 1: Simulation Data

Table 2 shows the variables and values used in the simulation of the well shown in Figure 5. The other variables used in the simulation are functions of the values shown in the table.

6.3. Simulation 1: General Results

Figure 6, on page 12, shows the time evolution of the fluid temperature along the pipe. It can be seen that the dynamic change in the curve slows down as the operating time increases. In particular, the curve for 2 hours differs little from the curve for 6 hours.

Figure 7, on page 13, shows the time evolution of the radial temperature at a depth of 5 meters. As can be seen in the graph, after 6 hours, the thermal disturbance was approximately restricted to a distance of 50 centimeters from the center of the pipe, validating the external radius of the reservoir used in the simulation (1 meter). Given the simulation time considered (6 hours), using a small value for the outer radius of the reservoir (Reserv_Radius) would result in overheating/undercooling of the well, saturating the simulation and making the data generated unrepresentative of a real situation.

Therefore, for the prototype installed at UFRN, which has an outer radius of approximately 30 centimeters, we recommend adopting 2 hours as the maximum operating time. After this time, the prototype will overheat, and the data generated will no longer be valid.

6.4. Simulation 1: Inferred Flow

The graph in Figure 8, on page 14, shows the evolution of the flow measurement for the methods and functions presented in this work. As can be seen, all the methods tend to converge to the reference flow used in the simulation (

0.3 * 10^{- 3} [m^{3} / s]

). However, the methods derived from the constant temperature assumption show better results, i.e., constant temperature at

r_{1}

(curve Constant Temp), constant fluid temperature with radiation condition at

r_{1}

(curve Radiation Condition) and the Hagoort approximations (curve Hagoort). This greater accuracy of the constant temperature methods is a consequence of the short transit time of the fluid in the pipe, equal to 1 minute, which means that the fluid temperature at the point of measurement has a small time variation compared to the time variation of the heat flow.

The graph in Figure 9, on page 14, shows the evolution of the flow measurement error. As expected, the choice of transient function will result in a greater or lesser intensity systematic error. Applying the proposed analytical solution results in a significantly lower systematic error for the entire simulated operating time.

6.5. Simulation 2: Well with Simplified Completion

This second simulation aims to verify the impact of transient heat transfer in the well completion on the inference of flow rate using the methods and formulations outlined in this article. To achieve this, consider the well completion represented in Figure 10. This completion technique, which involves only tubing and cementation, is widely utilized in real injector and producer wells.

As in the previous simulation, adiabatic boundaries were used at the edges of the simulated well, making it easy to detect saturation. A constant temperature condition was also used at the inlet of the injection column (tubing).

6.6. Simulation 2: Simulation Data

Table 3 shows the variables and values used in this topic. Following the same approach as the previous simulation, the other variables used are functions of the values presented.

6.7. Simulation 2: General Results

The graph in Figure 11, on page 17, shows the time evolution of the fluid temperature along the pipe. As in the previous simulation, it can be seen that the dynamics of the change in the curve lose speed as time progresses, but with less intensity, resulting in higher temperatures at the end of each measurement time. As observed in the previous simulation, the fluid temperature curve for 2 hours of operation differs little from the curve for 6 hours.

The graph in Figure 12, also on page 17, shows the time evolution of the radial temperature at a depth of 5 meters, half the length of the pipe. After 6 hours, the thermal disturbance was restricted to a distance of 50 centimeters from the center of the pipe, which shows that the simulation was not saturated for the time considered.

6.8. Simulation 2: Inferred Flow

When flow was measured using classical methods, simplifications were made in the calculation of the global heat transfer coefficient, i.e. the transient of heat transfer through the completion (tubing and cementing) was not considered. When measured using the proposed analytical method, only the transient effects of the tubing were ignored due to its high thermal conductivity when compared to the other layers of the well (cementing and formation).

The graph in Figure 13, on page 18, shows the evolution of flow measurement for the methods and transient functions presented in this article.

As can be seen in the Figure 13, as in the previous simulation, all the methods tended to converge to the reference flow rate (

0.3 * 10^{- 3} [m^{3} / s]

). However, unlike the previous simulation, the methods derived from the assumption of constant heat flow, i.e., constant heat flow at

r_{1}

(curve Constant Flow) and the Hasan & Kabir approximations (curve Hasan), showed better results. The method based on applying the radiation boundary condition in conjunction with the assumption of constant fluid temperature (Radiation Condition curve) showed good results in both simulations, given its flexibility in dealing with heat flow into the reservoir. In the same sense, the proposed analytical solution (curve Proposed Solution) gave flexibility to the fluid temperature and showed visibly better results than all the other methods studied.

The graph in Figure 14 shows the evolution of the systematic error in flow measurement. The systematic error resulting from applying the proposed analytical method (curve Proposed Solution) tends to be less than 1% in the first moments of the simulation, unlike the other methods, which require hours or days of operation. The systematic error in the first instants of the analytical approach is a consequence of neglecting the heat flow in the vertical direction, which is significant in the first instants in the vicinity of the well, and of not taking into account the effects of the pipe on heat propagation in the system as a whole.

7. Conclusions and Future Work

As we have seen, the quasi-static methods have a systematic error inherent to the methodology, which results from ignoring the variation of the internal energy of the fluid in the tubing and, above all, from ignoring the transience of the heat transfer in the completion. For the first few hours of well operation, this systematic error can be higher than 20%, depending on the completion and the transient function chosen. On the other hand, considering a simplified completion (tubing and cementing), the proposed analytical approach, which is the result of solving the heat transfer differential equations, shows a systematic error of <1% in the first two hours of operation.

In terms of methodology, the results of this study demonstrated a mathematical and computational way to improve the quality of flow measurement in injection wells in the early stages of operation. For this purpose, it was necessary to take into account the transient nature of the heat propagation in the completion as well as the temporal variation of the fluid temperature.

As a complement to the research carried out in this paper on reducing systematic error in flow measurement in injection wells, the following studies and actions can be developed in future work:

(a): to model and to solve the system of differential equations for the standard completion (tubing, annulus, casing and cementation);
(b): to model and to develop the computer code to include the effects of combined convection (natural convection and radiation) in the annulus;
(c): to study the effects of vertical heat flow in the completion, if we are looking for greater anticipation in the convergence of the calculated flow values;
(d): to research and to develop techniques to minimize the dependence on initial conditions.

Author Contributions

V. S. L., C. W. S. d. P. M. and A.O.S. conceived and designed the study; V. S. L., C. W. S. d. P. M. and G. A. E. E., were responsible for the methodology; V. S. L., G. P. d. O., D. A. d. M. F. and W. L. A. d. S. performed the simulations and experiments; V. S. L., E.R.L.V., G. P. d. O., D. A. d. M. F. and W. L. A. d. S. reviewed the manuscript and provided valuable suggestions; V. S. L., E.R.L.V., G. P. d. O., D. A. d. M. F. and W. L. A. d. S. wrote the paper; G. A. E. E., C. W. S. d. P. M., and A.O.S. were responsible for supervision. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded in part by Universidad Nacional de San Agustin de Arequipa (UNSA), and supported by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001 and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

$α$	Thermal diffusivity $[m^{2} / s]$
$α_{1}$	Thermal diffusivity of medium 1 $[m^{2} / s]$
$α_{2}$	Thermal diffusivity of medium 2 $[m^{2} / s]$
$α_{3}$	Thermal diffusivity of medium 3 $[m^{2} / s]$
$β_{1}$	Biot number relative to the medium 1
$β_{f o r}$	Relative Biot number of the formation
$μ$	Absolute viscosity $[N . s / m^{2}]$
$θ_{1}$	Dimensionless temperature of the medium 1
${\bar{θ}}_{1}$	Dimensionless temperature of medium 1 in the Laplace domain
$θ_{f}$	Dimensionless injection fluid temperature
${\bar{θ}}_{f}$	Dimensionless temperature of the injection fluid in the Laplace domain
$ρ$	Mass density $[k g / m^{3}]$
$ρ_{0}$	Mass density of the medium 0 $[k g / m^{3}]$
$ρ_{1}$	Mass density of the medium 1 $[k g / m^{3}]$
$ρ_{2}$	Mass density of the medium 2 $[k g / m^{3}]$
$ρ_{3}$	Mass density of the medium 3 $[k g / m^{3}]$
$ρ_{f o r}$	Mass density of the formation $[k g / m^{3}]$
$τ$	Shear stress $[N / m^{2}]$
$τ$	Dimensionless Fourier time
$a_{w}$	Angular coefficient of geothermal temperature with respect to w $[^{\circ} C]$
$b_{w}$	Linear coefficient of geothermal temperature with respect to w $[^{\circ} C]$
c	Specific heat $[J / (k g . K)]$
$c_{0}$	Specific heat of the medium 0 $[J / (k g . K)]$
$c_{1}$	Specific heat of the medium 1 $[J / (k g . K)]$
$c_{2}$	Specific heat of the medium 2 $[J / (k g . K)]$
$c_{3}$	Specific heat of the medium 3 $[J / (k g . K)]$
$c_{f o r}$	Specific heat of formation $[J / (k g . K)]$
f	Friction factor
$f (t)$	Transient formation function
h	Heat convection coefficient $[W / (m^{2} . K)]$
$h_{c o m b}$	Combined heat transfer coefficient $[W / (m^{2} . K)]$
$I_{0}$	Modified Bessel function of first type and order 0
k	Thermal conductivity $[W / (m . K)]$
$k_{0}$	Thermal conductivity of the medium 0 $[W / (m . K)]$
$k_{1}$	Thermal conductivity of the medium 1 $[W / (m . K)]$
$k_{2}$	Thermal conductivity of the medium 2 $[W / (m . K)]$
$k_{3}$	Thermal conductivity of the medium 3 $[W / (m . K)]$
$k_{c i m}$	Thermal conductivity of cementation $[W / (m . K)]$
$k_{f o r}$	Thermal conductivity of the formation $[W / (m . K)]$
$k_{t u b}$	Pipe thermal conductivity $[W / (m . K)]$
$k_{r e v}$	Casing thermal conductivity $[W / (m . K)]$
$K_{0}$	Modified Bessel function of second type and order 0
L	Length (depth) of the well $[m]$
$\dot{m}$	Mass flow $[k g / s]$
r	Radial coordinate $[m]$
$r_{1}$	Pipe inner radius $[m]$
$r_{2}$	Pipe external radius $[m]$
$r_{3}$	Casing inner radius $[m]$
$r_{4}$	Casing outer radius $[m]$
$r_{5}$	Cementation external radius $[m]$
$R e$	Dimensionless Reynolds number
$P r$	Prandtl dimensionless number
$\dot{Q}$	Heat transfer rate $[W]$
s	Transformation variable for the Laplace domain
t	Time $[s]$
T	Temperature $[^{\circ} C]$
$T_{f}$	Average fluid temperature $[^{\circ} C]$
$T_{g}$	Geothermal temperature $[^{\circ} C]$
U	Overall heat transfer coefficient $[W / (m^{2} . K)]$
$\tilde{U}$	Average overall heat transfer coefficient $[W / (m^{2} . K)]$
$v_{m e d}$	Average velocity of the fluid velocity profile $[m / s]$
w	Dimensionless vertical coordinate
z	Vertical coordinate $[m]$
$z_{1}$	Initial vertical coordinate of the region under analysis $[m]$
$z_{2}$	Final vertical coordinate of the region under analysis $[m]$
$\tilde{z}$	Average vertical coordinate $[m]$

References

Thomas, J. E. et al. Fundamentos de engenharia de petróleo. 2. ed. Rio de Janeiro: Interciência, 2001.(In Portuguese).
Reges, J. E. O.; Salazar, A. O.; Maitelli, C. W. S. P.; Carvalho, L. G.; Britto, U. J. B. Flow Rates Measurement and Uncertainty Analysis in Multiple-Zone Water-Injection Wells from Fluid Temperature Profiles. Sensors 2016, 16(07), 1077. [Google Scholar] [CrossRef] [PubMed]
Wu, Y. S.; Pruess, K. An analytical solution for wellbore heat transmission in layered formations. SPE Reservoir Engineering 1990, 05(04), 531–538. [Google Scholar] [CrossRef]
Ramey Jr, H.; Henry, J. Wellbore heat transmission. Journal of petroleum Technology 1962, 14(04), 427–435. [Google Scholar] [CrossRef]
Assmann, B. W. Previsão do comportamento de pressão e temperatura transitorios em poços de petroleo e oleodutos. Dissertação (Mestrado) - Universidade Estadual de Campinas. Campinas: UNICAMP, 1993. xii, 128 p.(In Portuguese).
Lima, V. S. Perfis térmicos em poços injetores d’água: redução do erro sistemático na medição de vazão durante o regime transitório. Ph.D. Thesis, UFRN University, Natal, Brazil, 2019. (In Portuguese). [Google Scholar]
Çengel, Y. A.; Ghajar, A. J. Tranferência de calor e massa: uma abordagem prática. 4. ed. Porto Alegre: AMGH, 2012. xxii, 904 p.(In Portuguese).
Hasan, A. R.; Kabir, C. S. Aspects of wellbore heat transfer during two-phase flow. SPE Production & Facilities 1994, 9, 211–216. [Google Scholar]
Matthews, C. S.; Russell, D. G. Pressure buildup and flow tests in wells. New York: Henry L. Doherty Memorial Fund of AIME New York, 1967. v. 1.
Hagoort, J. Ramey’s wellbore heat transmission revisited. SPE journal, Society of Petroleum Engineers 2004, 9, 465–474. [Google Scholar] [CrossRef]
Jacquot, R.; Steadman, J.; Rhodine, C. The gaver-stehfest algorithm for approximate inversion of laplace transforms. IEEE Circuits & Systems Magazine 1983, 5, 4–8. [Google Scholar]

Figure 1. Standard well schematic.

Figure 2. Control volumes and conservation of energy.

Figure 3. Algorithm for solving the equation system.

Figure 4. New Control volumes and conservation of energy.

Figure 5. Well without completion - details of the problem.

Figure 6. Well without completion - Longitudinal temperature evolution.

Figure 7. Well without completion - Radial temperature evolution.

Figure 8. Well without completion - Measurement at 9 meters.

Figure 9. Well without completed - Measurement error at 9 meters.

Figure 10. Completion with two layers - Details of the problem.

Figure 11. Completion with two layers - Longitudinal temperature evolution.

Figure 12. Completion with two layers - Radial temperature evolution.

Figure 13. Completion with two layers - Measured at 9 meters

Figure 14. Completion with two layers - Measurement error at 9 meters.

Table 1. System of differential equations - Simplified completion.

	Description	Equation
Formation( $r_{3} \leq r \leq \infty$ )	General Equation	$\frac{\partial^{2} T (r, z, t)}{\partial r^{2}} + \frac{1}{r} \frac{\partial T (r, z, t)}{\partial r} = \frac{1}{α_{3}} \frac{\partial T (r, z, t)}{\partial t}$ , $α_{3} = \frac{k_{3}}{ρ_{3} c_{3}}$
	Initial Condition	$T (r, z, 0) = T_{g} (z)$
	Boundary Condition (∞)	$\underset{r \to \infty}{l i m} T (r, z, t) = T_{g} (z)$
	Boundary Condition ( $r_{3}$ )	a) $k_{2} \frac{\partial T (r_{3}^{-}, z, t)}{\partial r} = k_{3} \frac{\partial T (r_{3}^{+}, z, t)}{\partial r}$
		b) $T (r_{3}^{-}, z, t) = T (r_{3}^{+}, z, t)$
Cementation( $r_{2} \leq r \leq r_{3}$ )	General Equation	$\frac{\partial^{2} T (r, z, t)}{\partial r^{2}} + \frac{1}{r} \frac{\partial T (r, z, t)}{\partial r} = \frac{1}{α_{2}} \frac{\partial T (r, z, t)}{\partial t}$ , $α_{2} = \frac{k_{2}}{ρ_{2} c_{2}}$
	Initial Condition	$T (r, z, 0) = T_{g} (z)$
	Boundary Condition ( $r_{3}$ )	a) $k_{2} \frac{\partial T (r_{3}^{-}, z, t)}{\partial r} = k_{3} \frac{\partial T (r_{3}^{+}, z, t)}{\partial r}$
		b) $T (r_{3}^{-}, z, t) = T (r_{3}^{+}, z, t)$
	Boundary Condition ( $r_{2}$ )	a) $k_{1} \frac{\partial T (r_{2}^{-}, z, t)}{\partial r} = k_{2} \frac{\partial T (r_{2}^{+}, z, t)}{\partial r}$
		b) $T (r_{2}^{-}, z, t) = T (r_{2}^{+}, z, t)$
Tubing( $r_{1} \leq r \leq r_{2}$ )	General Equation	$\frac{\partial^{2} T (r, z, t)}{\partial r^{2}} + \frac{1}{r} \frac{\partial T (r, z, t)}{\partial r} = \frac{1}{α_{1}} \frac{\partial T (r, z, t)}{\partial t}$ , $α_{1} = \frac{k_{1}}{ρ_{1} c_{1}}$
	Initial Condition	$T (r, z, 0) = T_{g} (z)$
	Boundary Condition ( $r_{2}$ )	a) $k_{1} \frac{\partial T (r_{2}^{-}, z, t)}{\partial r} = k_{2} \frac{\partial T (r_{2}^{+}, z, t)}{\partial r}$
		b) $T (r_{2}^{-}, z, t) = T (r_{2}^{+}, z, t)$
	Boundary Condition ( $r_{1}$ )	a) $- r_{1}^{+} \frac{\partial T (r_{1}^{+}, z, t)}{\partial r} = β_{1} [T_{f} (z, t) - T (r_{1}^{-}, z, t)]$ , $β_{1} = \frac{r_{1} h}{k_{1}}$
		b) $T (r_{1}^{-}, z, t) = T (r_{1}^{+}, z, t)$
Fluid( $0 \leq r \leq r_{1}$ )	General Equation	$\frac{\partial T_{f} (z, t)}{\partial t} + v_{m e d} \frac{\partial T_{f} (z, t)}{\partial z} = \frac{2 h}{r_{1} ρ_{0} c_{0}} [T (r_{1}, z, t) - T_{f} (z, t)]$
	Initial Condition	$T_{f} (z, 0) = T_{g} (z)$
	Boundary Condition ( $z_{0}$ )	$T_{f} (0, t) = T_{f 0}$
	Boundary Condition ( $r_{1}$ )	a) $- r_{1}^{+} \frac{\partial T (r_{1}^{+}, z, t)}{\partial r} = β_{1} [T_{f} (z, t) - T (r_{1}^{-}, z, t)]$ , $β_{1} = \frac{r_{1} h}{k_{1}}$
		b) $T (r_{1}^{-}, z, t) = T (r_{1}^{+}, z, t)$

Table 2. Well without completion - Data used in the simulation (L=10 meters).

	Variable	Value	Description
General data	dt0	0.25 s	Time step
	dz0	0.1 m	Space step in the z direction (longitudinal)
	tsimu	6 h	Simulation time
	p_samples_z	0.1 m	Sampling period in the "z" direction
	p_samples_t	15 s	Time sampling period
	TAmb	35 $^{\circ} C$	Environment temperature
	TgeoType	1	Geothermal temperature type (0=constant; 1=linear)
	ThermalSource	0	0=Adiabatic boundary; 1=Thermal source equal to Tgeo
Geometric	L0	10 m	Longitudinal length of the well to be simulated
	Tub_Radius	0.0254 m	[r1] Pipe inner radius
	Reserv_Radius	1 m	[r2] Reservoir external radius
	DIV0	1	Number of region radial divisions 0
	DIV1	40	Number of region radial divisions 1
Fluid	k0	0.636 $W / (m . K)$	Water thermal conductivity
	Cp0	4184 $J / (k g . K)$	Water-specific thermal capacity
	ro0	1000 $k g / m^{3}$	Water-specific mass
	mi0	0.0006 $N . s / m^{2}$	Absolute viscosity
	f0_	0.0003 $m^{3} / s$	Volume flow
	TF_IN	20 $^{\circ} C$	Inlet fluid temperature
Formation	k1	2.42 $W / (m . K)$	Reservoir thermal conductivity
	Cp1	1500 $J / (k g . K$ )	Reservoir-specific thermal capacity
	ro1	2100 $k g / m^{3}$	Reservoir specific mass
	GradGeo	0.365 $^{\circ} C / m$	Geothermal gradient in the "z" dimension
	TSurf	TAmb	Surface temperature

Table 3. Completion with two layers - Data used in the simulation.

	Variable	Value	Description
General data	dt0	0.25 s	Time step
	dz0	0.1 m	Space step in the z direction (longitudinal)
	tsimu	6 h	Simulation time
	p_samples_z	0.1 m	Sampling period in the "z" direction
	p_samples_t	15 s	Time sampling period
	TAmb	35 $^{\circ} C$	Environment Temperature
	TgeoType	1	Geothermal temperature type (0=constant; 1=linear)
	ThermalSource	0	0=Adiabatic boundary; 1=Thermal source equal to Tgeo
Geometric	L0	10 m	Longitudinal length of the well to be simulated
	Tub_Radius	0.0254 m	[r1] Longitudinal length of the well to be simulated
	Reserv_Radius	1 m	[r4] Reservoir external radius
	DIV0	1	Number of region radial divisions 0
	DIV1	3	Number of region radial divisions 1
	DIV2	10	Number of region radial divisions 2
	DIV3	30	Number of region radial divisions 3
Fluid	k0	0.636 $W / (m . K)$	Water thermal conductivity
	Cp0	4184 $J / (k g . K)$	Water-specific heat capacity
	ro0	1000 $k g / m^{3}$	Water specific mass
	mi0	0.0006 $N . s / m^{2}$	Absolute viscosity
	f0_	0.0003 $m^{3} / s$	Volume flow
	TF_IN	20 $^{\circ} C$	Inlet fluid temperature
Tub	k1	14 $W / (m . K)$	Reservoir thermal conductivity
	Cp1	502 $J / (k g . K)$	Reservoir-specific thermal capacity
	ro1	8000 $k g / m^{3}$	Reservoir specific mass
	Tub_Thickness	0.635 $c m$	Pipe thickness
Cement	k2	0.9 $W / (m . K)$	Reservoir thermal conductivity
	Cp2	900 $J / (k g . K)$	Reservoir-specific thermal capacity
	ro2	2400 $k g / m^{3}$	Reservoir -specific mass
	Cim_Thickness	5.08 $c m$	cementation thickness
Formation	k3	2.42 $W / (m . K)$	Reservoir thermal capacity
	Cp3	1500 $J / (k g . K)$	Reservoir-specific thermal capacity
	ro3	2100 $k g / m^{3}$	Reservoir -specific mass
	GradGeo	0.365 $^{\circ} C / m$	Geothermal gradient in the "z" dimension
	TSurf	TAmb	Surface temperature

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Thermal Profiles in Water InjectionWells: Reduction of the Systematic Error in the Flow Measurement during the Transient Regime

Abstract

1. Introduction

2. Preliminaries

3. Quasi-Static Approach

4. Analitical Approach

5. Computer Simulation

6. Results

6.1. Simulation 1: Well without Completion

6.2. Simulation 1: Simulation Data

6.3. Simulation 1: General Results

6.4. Simulation 1: Inferred Flow

6.5. Simulation 2: Well with Simplified Completion

6.6. Simulation 2: Simulation Data

6.7. Simulation 2: General Results

6.8. Simulation 2: Inferred Flow

7. Conclusions and Future Work

Author Contributions

Funding

Conflicts of Interest

Abbreviations

References

MDPI Initiatives

Important Links

Subscribe