Laplace transform-based modeling of heat transfer in producing wellbores

This paper presents a novel approach to modeling heat transfer in hydrocarbon-producing wells, aiming to simulate annular pressure buildup (APB) in offshore wells. The proposed framework enhances traditional models by incorporating previously overlooked terms in pseudo-steady-state simulators. A Laplace transformation-based formulation of the governing differential equations is developed and solved to achieve better simulation results for shorter timespans. Two actual oil-producing wellbores are used for model validation. The first wellbore, a vertical well with the first annulus partially filled with N2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\hbox {N}_{2}}$$\end{document}, assesses the new formulation’s performance compared to the traditional pseudo-steady-state formulation. Additionally, a convergence test is conducted in this well to verify the level of additional computational complexity associated with the new formulation. Simulation results show that slower annuli heating leads to a slower APB, resulting in an 8 MPa difference in predicted pressurization. The second case, from a well undergoing an extended well test, provides field data to simulate wellbore heating due to heat transfer from the production stream. A comparison of the new formulation with the field data demonstrates improved temperature representation within the first few days of production, with an average deviation of 4 K in predicted wellhead-produced fluid temperature.


Introduction
Accurately simulating the temperature profile during hydrocarbon production operations in wells is essential for designing the wellbore structure correctly.As production begins and hot reservoir fluids flow through the production string, heat is exchanged with the surrounding rock formation due to a geothermal gradient of approximately 25 Kelvin per kilometer drilled in the earth's crust [1].Because drilling involves creating annular spaces between casings and adjacent tubing, every wellbore contains trapped fluids.As heat flows radially outward from the production stream, the annular fluids heat up, expand, and experience a pressure increase.Failure to account for the annular pressure buildup (APB) during wellbore design can lead to structural collapse and a premature end to the well's life cycle [2].
APB has garnered significant attention in the research community ever since its correlation with incidents in the field was confirmed.One notable instance, documented by Bradford et al. [3], involved a production casing collapsing under excessive load shortly after the commencement of well production.Other researchers, such as [4][5][6][7], have also delved into investigating this failure.
To ensure an accurate prediction of pressure increase resulting from APB, it is essential to have precise estimations of heat transfer.The seminal paper of Ramey [8] established the fundamental hypotheses and governing equations for heat transfer in wellbores, which remain the basis for more advanced simulators [9].In general, the formation is considered a semi-infinite medium, and the impact of wellbore thermal capacitance becomes insignificant after a few days of production.Moreover, due to the wellbore's slenderness (length-to-diameter ratio), onedimensional semi-steady-state approaches are commonly Technical Editor: Francis HR Franca.
652 Page 2 of 19 used to solve the momentum and energy balances, where the heat flow term varies only in the radial direction.
Over the years, the solution procedure and initial hypotheses proposed by Ramey [8] have been refined and enhanced by other authors.For instance, one notable improvement was made by Durrant and Thambynayagam [10], who were among the first to introduce the marching integration procedure for the wellbore energy equation.This advancement enabled the consideration of temperature and pressure dependence of the annular fluid properties.Additionally, the marching scheme facilitated the formulation of the radial heat flow term, particularly in complex multi-string wellbore geometries.
Several authors have examined different aspects of the problem.Some focused on the occurrence of multiphase flow in the production tubing [11][12][13][14], while others investigated the various heat transfer mechanisms in the concentric annuli [15].Furthermore, several methods have been developed to solve the heat diffusion equation in the surrounding formation [9,[16][17][18][19].In particular, Hagoort [9] addressed the limitations of the pseudo-steady-state approach, which tends to overpredict the wellbore thermal profile during short-term operations.
Regarding APB prediction, Halal and Mitchell [20] proposed a coupled formulation that relies on the deformation of the multi-string casing.Additionally, Oudeman and Bacarreza [21] and Oudeman and Kerem [22] investigated the mechanisms of pressure buildup resulting from thermally driven fluid expansion.They introduced a general linearized method to quantify the APB.
Recently, da Veiga et al. [23] and Martins et al. [24] proposed different thermo-structural wellbore models that incorporate multiple physics.The latter model, which is a transient two-dimensional approach, considers thermal capacitance and longitudinal heat conduction in the wellbore components and surrounding formation.This makes it particularly well-suited for short-term operations like liquid or gas injection.On the other hand, the former model is a computationally inexpensive one-dimensional pseudo-steady approach, making it suitable for real-time simulations.However, using the pseudo-steady model without considering wellbore thermal capacitance effects limits its application to situations with low Fourier numbers (early time).
The present paper aims to bridge the gap between two formulations by introducing a novel method that enhances the prediction of wellbore thermal behavior throughout all stages of oil production.In contrast to traditional approaches [23], which exclude thermal capacitance terms and focus on pseudo-steady characteristics for rapid computation of temperature and pressure fields, our proposed radial heat transfer model incorporates these previously neglected terms.To solve the model, we utilize a hybrid Laplace transform method and validate its accuracy using field data from two real wellbores.
The first wellbore chosen for validation is the same as the one presented in [23].This particular wellbore is involved in oil production and features a partially filled first annulus with N 2 .The availability of data from previous operations makes it crucial for assessing the efficacy of our new formulation in handling long-term operations, which can be accurately represented using pseudo-steady-state simulators.
For the second wellbore, we utilize data from the initial production phase to evaluate the new formulation's ability to better represent the thermal behavior of the produced fluid during shorter operations.Specifically, we focus on the temperature of the produced fluid within the first couple of days of production.

Problem formulation
The general structure of a well completed with three annuli is depicted in Fig. 1.Two of the annuli are filled with drilling mud or completion fluid, while the third annulus is cemented up to the wellhead.The thermal problem involves solving differential equations along the well.These equations can be divided into two parts: one concerning the production fluid, which includes two-phase momentum and energy balances, and the other concerning the sealed annuli, which involves solving equations for mass conservation and hydrostatic pressure relationships.
The conservation equations for the production stream assume a pseudo-steady-state form by neglecting the transient and kinetic energy terms in the momentum and energy balances, as well as the heat conduction in the axial direction in the energy equation [23].The equations can be expressed as follows: In these equations, P represents the pressure, f s denotes the two-phase friction factor (considering slip), ns is non-slip fluid density at the local P and T, v m signifies the mixture velocity, and d T corresponds to the inner diameter of the pro- duction tubing.The values of f s and s are calculated using the Beggs and Brill [25] correlations, incorporating the corrections proposed by Brill and Beggs [26].Furthermore, h represents the enthalpy of the produced fluid, ṁ denotes the mass flow rate, g signifies the acceleration due to gravity, represents the inclination from the horizontal, and Q ′ refers to the heat flow per unit length.In both equations, z is the wellbore length.As most wells deviate from the vertical, the integration path follows the curvature of the well.
In each sealed annulus, the local pressure is calculated by solving the hydrostatic pressure and mass conservation equations as follows: where r out and r inn are the local outer and inner annulus radii.By writing the mass conservation in this format, it is possible to couple the structural deformation models with the thermally driven fluid expansion effects by accounting for the casing radial displacements.
In order to determine the local density of the packer fluids, it is necessary to have knowledge of the local pressure and temperature.While the former can be obtained by integrating Eq. ( 3), the latter is calculated using an equivalent thermal network that determines the heat flow per unit length, as described in Eq. (2).Consequently, at each depth, a corresponding thermal network is constructed by taking into account the thermal capacitances of the materials involved, such as casings, trapped fluids, and cement.This assembly is depicted in Fig. 2, which illustrates the thermal network for the specific region highlighted in Fig. 1.
In Fig. 2, the subscripts indicate which interface was used to calculate the temperatures and heat fluxes.For instance, T 0 represents the bulk temperature of the production stream, while T 1 is the temperature at the wall of the production tub- ing.In the geometry shown in Fig. 1, the production tubing is encased by an annulus filled with packer fluid (between T 2 and T 3 ), followed by another steel tube (between T 3 and T 4 ) and finally a cement layer in contact with the rock formation.
At each interface between adjoining media, an energy balance accounts for the thermal energy stored in the material surrounding the preceding node.Thus: where the heat rate per unit length at each node is calculated based on the thermal resistance per unit length as follows: The resistances per unit length are calculated based on heat transfer coefficients, which rely on the heat transfer mode occurring at the specific interface under analysis.For instance, the production fluid utilizes the forced convection correlation from Chen [27], excluding the contribution of nucleate boiling.When it comes to the fluids trapped in the annuli, the thermal resistance consists of two components in parallel.One component is the resistance caused by natural convection, which is calculated using the correlation Fig. 2 Equivalent resistances and capacitances of the thermal network (red contour region of Fig. 1.).Here, the index 0 is equivalent to the produced fluid, and the index " N m " is the last node representing the contact interface from the wellbore and the rock formation, also called "form, int" 652 Page 4 of 19 provided by Zhou [28].The other component is the resistance resulting from radiation between the solid surfaces, as suggested by Hasan and Kabir [15].Lastly, the cylindrical walls, namely the cement, tubing, and casings, are governed by expressions for pure radial conduction resistances [29].
The heat conduction at the outermost node, where the wellbore meets the formation, is solved as a two-dimensional boundary value problem (BVP).This is coupled with the thermal network through the following equation: where r fm is the formation inner radius, k fm is the formation thermal conductivity, and T fm is the formation temperature.Q ′ fm is the heat flow per unit length at the outermost node of Fig. 2, i.e., where T N m = T form,int .
The governing equation for the formation temperature in the BVP assumes an isotropic semi-infinite medium and can be expressed as follows: with the following boundary conditions where T int represents the temperature at the outermost node of the thermal network, and T init is the initial temperature of the boundary value problem, which in this case is the geothermal temperature.fm denotes the thermal diffusivity of the rock formation.
To non-dimensionalize the governing equations, the following dimensionless groups are introduced: and The dimensionless BVP is resolved using the Laplace transform technique, wherein the resulting algebraic system is addressed using simpler methods [30].The solution in the Laplace domain can be expressed as follows: where r D represents the dimensionless radius, T D corre- sponds to the dimensionless temperature in the Laplace domain, and s is the time variable transformed into the Laplace domain.The function K 0 refers to the modified Bes- sel function of the second kind with order 0. By applying the Laplace transform to each equation representing the nodes in Fig. 2, it is possible to rewrite the set of differential equations as a simple linear system to be solved in the Laplace domain.In the first node, i = 0 , the bulk temperature of the production fluid is prescribed.Note that this temperature is obtained from the solution of the conservation of energy in the produced fluid.Thus: The energy balance at the interface between the production stream and the tubing, i = 1 , neglects the thermal capaci- tance of the fluid due to a comparatively higher magnitude of energy advection.Thus, without the transient term in Eq. ( 5), the equation for the temperature at the wall of the production tubing becomes: For a generic node between the second and the last nodes, i ∈ {2, 3, … , N m − 1} , the general energy balance can be applied: Finally, the energy balance at the interface between the wellbore and formation is given by: Equations ( 15)-( 18) can be solved by simply inverting the coefficient matrix for a range of values of s in the Laplace domain.To convert the solution back to the time domain, numerical inversion techniques can be employed.In this paper, the concentrated matrix exponential (CME) method [31] was utilized.This method is based on the Abate-Whitt family of numerical inversion methods [32].Although several methods exist within this framework, the CME offers greater versatility in the inversion step.For instance, the algorithm proposed by Stehfest [33] is usually faster but limited to a small number of inversion nodes.On the other hand, ( 14) the CME provides hundreds of inversion nodes, thereby increasing the likelihood of finding a converged solution.
The subsequent stage of the numerical solution includes calculating thermo-elastic deformations related to each casing, tubing, cement, and the formation itself.In this study, we utilized the identical formulation proposed by da Veiga et al [23], which is based on the model introduced by Halal and Mitchell [20].This formulation employs an axisymmetric approach to depict the plane deformation at each deformation node.The key benefit of utilizing this model is its simplicity in integrating with the procedure and maintaining a consistent formulation that ensures stress and deformation continuity, even when multiple solids are in contact.
Equations ( 15)-( 18) are solved at each point along the integration of the momentum, energy and mass balance equations, Eqs. ( 1)-( 4), in the z-direction.After solving the entire wellbore, the total mass of each annulus is calculated, allowing for a comparison with its initial mass.The difference in the results is utilized to establish a residue equation, which can be iteratively solved to determine the correct pressurization of each annulus.The procedure involves modifying the initial pressure condition in each annulus until the mass matches the system's initial mass after heating.A multidimensional root finder is employed to solve this problem.
The boundary conditions for Eqs. ( 1) and ( 2) are determined using data on the pressure and temperature of the production fluid at the Permanent Downhole Gauge (PDG).For each annulus, the boundary conditions require the mass at the bottom to be zero, and the equations are then integrated to calculate the total mass based on the local fluid density and radial dimensions.

The pseudo-steady-state formulation
It is crucial to verify that the Laplace-based approach accurately reproduces the traditional pseudo-steady-state formulation over sufficiently long periods.The first step involves examining the behavior of Eq. ( 5) in the longterm limit, which results in slight variations in the temperature time derivative.Once this occurs, the heat balance at each interface can be determined solely by the thermal network's resistance terms.The heat balance at any interface can be mathematically expressed as: This balance is applied to every node except the one in contact with the formation.However, because the temperature time derivative in the last node is approximately zero, the heat equation in the formation follows the same behavior as steady heat flow.
The solution for the heat flow entering the formation under these conditions is given by: Since this solution cannot be analytically inverted, many authors suggest approximation functions to estimate the term.These approximation functions are presented in several works [9,15,34].
In this paper, to compare the new model and the traditional pseudo-steady-state approach, we used the time function proposed by Hasan and Kabir [18].According to the literature, this is a reliable approximation for long-term operations.
Regarding the challenges associated with implementing both models, a comprehensive program has already been developed in a previous paper [23].This program employs the pseudo-steady-state formulation and undergoes incremental improvements, allowing for the choice of either formulation to solve heat transfer problems.The sole distinction in the implemented code lies in the section responsible for solving the radial thermal system.Instead of solving the thermal system using steady-state equations and the formation time function, the new system is dependent on the Laplace domain's time variable, and it necessitates multiple solutions for the thermal system based on the number of inversion nodes used.

Results
This section presents the results of simulations carried out on real wellbores.Firstly, the results from the vertical well, which was previously used in [23], will be presented.After that, the results from a real well in the extended well test (EWT) will also be discussed.

Vertical wellbore
This well is located at a sea depth of 1394 m and has a single production zone situated between the depths of 4390 m to 4564 m.Since the well has a nearly vertical orientation, with negligible deviations based on field data, there is no distinction made between True Vertical Depth (TVD) and Measured Depth (MD), and only the latter is presented in the figures.Figure 3 depicts the wellbore geometry and the surrounding lithology.
It is important to mention some characteristics of this wellbore, which add a little more complexity to the present study.Firstly, annulus A (the innermost annulus) is filled (20)  with two completion fluids: nitrogen and brine.This is due to the gas lift operations conducted in the well.The geothermal gradient of the wellbore is depicted in Fig. 4. The deviation from the linear trend, observed beyond 4000 m, can be attributed to the presence of an Anhydrite layer.Anhydrite is an evaporitic rock known for its high thermal diffusivity, which results in a steeper gradient in its vicinity.The thermophysical properties of the formation discussed in this paper have been compiled from various sources in the literature [35][36][37].
According to the recommendations by [35], the drilling fluid is assumed to have properties similar to a 50 wt.%glycerin/water solution.The brine is modeled as an aqueous solution of NaCl, with a mass fraction that matches the density under standard conditions [38].
A PVT flash analysis yielded the composition of the produced fluid, with knowledge of the oil's composition up to C 20 .This composition is then inputted into commercial software, specifically Multiflash v.7.2 [39], to generate thermodynamic tables required for property calculations.
The input data for the model were obtained through various methods.Pressure and temperature data obtained from sensors placed inside the well-provided instantaneous values throughout the production period.On the other hand, flow data are only available as time-averaged values, especially for the specific time interval considered in the study.Thus, a single constant value will be assumed for the flow.
To validate the models, it is necessary to address certain limitations imposed by the production history.Figure 5 illustrates the data obtained at the Temperature-Pressure Transducer (TPT) during production.It is evident that this well undergoes multiple short production intervals, which complicates the analysis using the pseudo-steady-state model.As a result, the highlighted red region, as described in [35], was selected for the previous studies.However, from a thermal modeling perspective, this is not ideal because the well has been in production for over 90 days, and the formation is likely already heated, potentially influencing the final results.This aspect must be considered when comparing the models.
Other essential sensors are positioned within the PDG, which is situated at a depth of 3958.15 m, to monitor the pressure and temperature within the production tubing.However, the data from these sensors are limited to a specific timeframe, indicating measurement failures at other points in time.Consequently, Fig. 6 displays the data obtained from the PDG, with the highlighted region representing the 10-day interval utilized for the simulations.
Figure 6 supports the notion that all production intervals preceding the 10-day window (red zone) continue to impact the new "start" of the well.This is supported by the fact that the temperature profile has not yet reached a stable state at the geothermal gradient following a period of approximately 10 days without production ( t = −10 to 0 days).Therefore, it is reasonable to assume that, upon resumption of the flow, the production fluid will pass through regions that are still warm.As a result, simulations with shorter production times should be regarded as approximations and are likely to exhibit a slower response compared to the field data.The region of study for the thermal model, as well as the input data at the PDG, is depicted in Fig. 7.The time intervals between measurements are set at 6 h, starting from the moment when the first pressure fluctuation was recorded in the PDG meter, which also marks the beginning of the production process.This interval lasts for 10 days, covering the entire period of production, including the transition to pseudo-steady-state in the annuli.The symbols in the figure correspond to the input data points used in the simulation, and they were obtained by interpolating the data collected in the field.The hydrodynamic model assumes a constant flow rate of 1501 m 3 d −1 , which represents the average for the production region, and a gas-liquid ratio (GLR) of 211 m 3 m −3 .To model the multiphase flow within the produc- tion column, the correlation developed by [25] was used.For comparison purposes, the time function proposed by [18] was employed to model wellbore heat transfer.
Before presenting the results, a brief convergence test is conducted for the model.The results of the convergence test are displayed in Fig. 8.In this figure, the abscissa represents the number of inversion nodes used in the inverse transform to the thermal network.On the ordinate, the convergence criterion is presented, which is based on the absolute relative deviation from the APB predicted by the previous number of nodes.Starting with 4 inversion nodes, the number of nodes is doubled at each iteration, and the value of APB in annulus A is checked.As a general rule, the method was considered converged when the deviation caused by doubling the number of nodes was smaller than 1%.This condition was met with 16 nodes.For the converged case, a threefold increase in computational cost was observed when running simulations on two separate computers.In the pseudo-steady-state simulation, each simulated time step took an average of 5 s to perform the mass convergence iteration, while the hybrid Laplace model required approximately 15 s.Given that these formulations are primarily intended for performing simulations at specific time points without the need for intermediate steps, the increase in computational cost is not significant for engineering applications.
The primary result generated by the models is presented in Fig. 9.It showcases the APB behavior in each annulus starting from the resumption of production until the tenth day of production.Several noteworthy observations can be made from this figure.Firstly, both models predict a low APB value for annulus A. This behavior can be attributed to the large amount of N 2 in annulus A, which results in reduced pressure variations compared to other annuli, owing to its greater compressibility.
Based on the findings presented in Fig. 9, it is evident that the rate of the APB calculated by the hybrid Laplace model is slower compared to the conventional approach.Moreover, the difference becomes more pronounced as the distance between the annulus and the produced fluid increases, i.e., the variation in APB is most significant in annulus C, reaching up to 12.86 MPa after six hours.
This behavior is expected and represents the main advantage of the hybrid Laplace model over the pseudo-steadystate model.The reason is that the pseudo-steady-state model predicts an almost instantaneous pressure increase at the beginning of production, without accounting for the necessary heating phase that must occur.In scenarios with time-dependent phenomena, such as the coupled modeling of creep-induced deformation [36], the pseudo-steady-state model's overprediction of APB can result in inaccurate estimations of APB in long-term operations.
The comparison between field data obtained from the TPT gauge and model predictions is presented in Fig. 10.The figure shows two distinct results.Firstly, both models were able to accurately reproduce the real data within a 2 K error for long enough times.Secondly, when analyzing the predicted pressure behavior in the wellhead, a clear offset of approximately 1 MPa can be observed.
There could be two primary factors that might have caused this deviation.The first reason, which is relatively straightforward, pertains to issues with measuring the oil flow rate.The basis for this suspicion arises from the insufficient information provided by the well operator regarding the location and method of measurement employed.The second factor could be associated with the two-phase flow correlations, specifically the application of [25] in this particular instance.One important factor that is evident from the data is the deviation of the hybrid Laplace model compared to the pseudo-steady-state approach at the beginning of production.It is worth noting in Fig. 6 that the temperature at the bottom hole is not stable before the production restarts.Consequently, it can be inferred that the temperature along the remaining wellbore section has not yet returned to the geothermal gradient.This has a consequence: The temperature detected by the TPT after production restarts reveals that the system still retains thermal energy from previous production operations, thereby compromising the data for the initial days of production.
To comprehend the effects of transient terms in heat transfer, the temperature profiles for each annulus are presented in Figs.11 and 12.The first image illustrates the profiles after 6 h of production, while the second image displays the temperature profiles after 10 days of production.
The temperature profile of each annulus after 6 h reveals two main behaviors.Firstly, there are discontinuities at each depth where new annuli are added or where there is a change in the wellbore material.Since both models are essentially 1-D in the radial direction, the introduction of new materials in the layer leads to step changes in the thermal resistances.Consequently, this causes discontinuities in the thermal profile.This behavior can only be accounted for (and corrected) if conduction in the axial direction is included in the formulation.
The second observation from Fig. 11 is the tendency for heating to occur in both the axial direction (from the bottomhole to the wellhead) and the radial direction (with the innermost annulus heating first before progressing outward to the other annuli).For instance, in annulus A, at a depth close to 2000 m, at the bottom of annulus C, a noticeable deviation is observed between both models.However, this deviation is less pronounced in the production tube and becomes even less significant as the MD increases.
To further understand these issues, Table 1 presents the average temperature difference between the pseudosteady-state model and the hybrid Laplace model.Positive values, thus, indicate that the pseudo-steady-state formulation yields higher average temperature values.The results confirm the previous statement regarding the heating process, which starts from the innermost annulus and progresses outward, as well as from the bottom upwards.This is evident from the smaller temperature differences observed in annulus A during all evaluated time intervals.
To assess the structural integrity of the wellbore, the American Petroleum Institute (API) recommendation was followed, stating that a maximum deformation of 5% is admissible for any tubular within the wellbore [40].Figure 13 presents the radial strain at the inner and outer solid walls of each annulus, allowing one to examine the influence of the thermal model on the wellbore structure.It can be observed that the strain does not pose a significant risk to the structure in both models.
When considering the structural formulation that incorporates thermo-elastic deformation, notable deformations occur at regions with discontinuous heat flux.Consequently, the most substantial deviation is observed at the top of the annuli, where the temperature difference is at its maximum.In both models, the highest strain is observed in annulus C, with values exceeding 2 × 10 −3 in the radial direction for the pseudo-steady-state formulation.
As depicted in Fig. 12, the temperature distribution is similar for both formulations after 10 days of production, resulting in converging strains.This convergence is further validated in Fig. 14, which presents the strain at the wellbore during the first 10 days of production.As expected, the strain predicted by both models is identical at every depth after 10 days.Additionally, even after 10 days, when the APB reaches its maximum values in the simulation and the annuli heating becomes stable, the wellbore deformation does not exceed 5% in any of the casings.Once again, the largest deformation, exceeding 3 × 10 −3 , is observed in annulus C.

Well in the extended well test
The geometry of the EWT well is depicted in Fig. 15, illustrating the MD and radial dimensions of the axisymmetric wellbore components and rock formation.The gray-scale regions represent the solid components, namely cement, casings, and the production tubing, while the blue-scale regions represent the fluids trapped in each of the three annuli.The innermost annulus, Annulus A, is filled with a 9.1-ppg brine.However, there is no information available regarding the composition of the drilling fluid used for the other annuli, except for their standard weight.As a representation of the trapped fluid, NaCl brine is used for both Annulus B and C, with respective densities of 9.6 ppg and 9.8 ppg.The EWT well was drilled in a predominantly shale rock formation at a sea depth of 2051 m.The production zone is located at a measured depth of 5857 m.Temperature and pressures are measured at PDG positioned at a depth of 5490.2 m and TPT at 2501 m .An approximately constant geothermal gradient is present in this wellbore.The temperature is 277.15K at the sea bed (MD = 2501 m) and 383.15 The field data for the EWT well represent the initial phase of its production life, allowing for a thorough exploration of the benefits of the new thermal formulation.This formulation takes into account the thermal capacitance effects in the wellbore structure, maximizing its capacity. Figure 16 illustrates the model's input parameters, including pressure and temperature data from the PDG, as well as the flow rate over time.The figure highlights two production intervals, with the second interval exhibiting a higher average production flow rate.By cross-correlating the two signals, a time lag of approximately 28 min was observed between the produced fluid flow rate and the PDG pressure data.The actual flow rate data are used as input for the model, requiring interpolation to match the acquisition frequency of the input data with the temporal resolution of the numerical model.The same procedure is applied to the TPT signals for pressure and temperature.
Severe oscillations in the flow rate data during the first twelve hours of production prevented the use of the raw data as a model input due to numerical instabilities in the mass conservation routine.To circumvent this problem, a constant flow rate of 615 m 3 d −1 was adopted during this initial period to ensure consistent behavior of the dependent variables.
For the first 6 h of production, a time step of 15 min was used in the numerical simulation.Subsequently, from that point until the end of the first production interval (15 days), the time step was increased to 6 h.The period between the two production intervals was not simulated.However, in the second production interval, specifically between days 18 and 29, a 12-hour time step was used.Therefore, a total of 104 individual simulations were performed using each wellbore heat transfer model (pseudo-steady-state and hybrid Laplace).
The wellhead temperatures and pressures predicted by each wellbore heat transfer model are presented in Fig. 17 and compared with the TPT field data.Both methods yielded satisfactory results when compared to the TPT data.During the initial 12 h of production, before the completion of the first 15-day period, the performances of the two approaches are similar.However, after the first day, the hybrid Laplace method outperforms the pseudo-steady-state method by a small margin.The differences between the results obtained from the models are less than 1 K, which falls below the measurement uncertainty for the field data.
During the second production interval, following a two-day cessation of production from the well, the hybrid Laplace method exhibited satisfactory performance.More importantly, the temperature behavior indicates that the hybrid approach, utilizing an average flow rate, remains relatively accurate even after a brief production interruption.However, employing the average flow rate approach in the second interval resulted in a less accurate representation of the pressure profile, as evidenced by an almost constant pressure deviation of 830 kPa.
Since the most significant difference in the predicted data occurs within the initial hours, Fig. 18 presents the same results, focusing on the first 12 h of production.Considering the challenges associated with obtaining stable flow rate readings within this time frame, the numerical results exhibit promising general trends, capturing some of the field data characteristics.Initially, both predicted temperatures are higher than the field data, with a difference of approximately 10 K for the hybrid Laplace model and 30 K for the pseudosteady-state model.However, as production progresses, both models tend to estimate a slower growth rate for the temperature.Consequently, at approximately 2 h of production, the hybrid Laplace model begins to under-predict the temperature, while the pseudo-steady-state model never reaches under-prediction but deviates slightly from the data after 4 h.As time elapses, the discrepancy between the models diminishes, as the influence of the initial thermal inertia of the well, which is accounted for in the hybrid Laplace model, fades.
During the first 3 h of production, both models accurately predicted the pressure behavior at the wellhead.At the 2-hour mark, the hybrid Laplace model achieved a prediction deviation of less than 10.95 kPa.Between 2 and 6 h, the hybrid Laplace model exhibited a maximum temperature deviation of 2 K in its predictions.In contrast, the pseudosteady-state model had an average deviation of 3.68 K, with a maximum deviation of 4.6 K.This discrepancy is particularly pronounced in regions where there is suspicion of inaccurate data measurement for oil flowrate, as both models fail to capture the rate of temperature growth, which is primarily controlled by the flowrate.To address this issue and improve temperature predictions at shorter times, it may be beneficial to explore alternative techniques for representing the flowrate in these regions where a lack of data appears to significantly impact the results.
In addition to the direct analysis of the field data, a simple statistical analysis can be employed to assess the deviation.Table 2 presents the average absolute deviation (AAD) of these models with respect to the TPT data.It is evident that, in general, both models exhibit a close-to-expected deviation of approximately 1.3% for the measured pressure.However, there is a more noticeable disparity in the results for temperature.The pseudo-steady-state model demonstrates an expected average deviation of 2.2%, whereas the hybrid Laplace model achieves 0.9%.Although this temperature difference may appear small, it becomes significant for Fig. 17 Comparison of the TPT temperature and pressure with predicted results from proposed models in the EWT well Fig. 18 Comparison of TPT temperature and pressure measurements in the EWT well with results from the two models for the first 12 h of production temperatures above 300 K, potentially leading to a deviation of over 4 K in the predicted wellhead temperature.Despite the seemingly minor variation, such disparities in wellhead temperature can result in larger deviations in APB due to the low compressibility of the trapped liquid.Additionally, the 4% difference is related to the average of all 104 data points, which includes simulations for extended operation times when the traditional approach exhibits only a small deviation from the field data.
To gain a deeper insight into the distribution of uncertainties among the 104 data points, the AAD was computed for both the first half-day of production and the period extending up to three days of production.During the initial 12 h, a total of 24 simulations were conducted for each method.In this scenario, the deviation from the field data for the pseudo-steady-state method was 7.41% for the temperature of the produced fluid (21.97 K).In contrast, the Laplace model achieved an AAD of 2.64% (7.91 K), indicating a superior representation of the initial temperature profile during production.
In the simulated time spanning up to three days of production, analysis was performed on 35 data points.In this context, the AAD for the pseudo-steady-state method was 5.54% (16.68 K), while for the Laplace model, it was 2.23% (6.93 K).In both cases, there was a reduction in the deviation from the wellhead temperature measured at the wellbore.However, the Laplace formulation proved to be more accurate when compared to the data.As previously mentioned, this behavior is expected since the Laplace formulation accounts for thermal capacitances, resulting in a slower heating of the wellbore.Consequently, it provides a more accurate representation of the initial heating during production.With the progression of production, both models should converge toward similar results, thus minimizing the disparities in the outcomes.
Figure 19 presents the APB at each annulus for the entire simulated interval.As can be seen, the pressure increase in Annulus A was the highest among all annuli of this well.The predicted pressure increase reached its maximum value at the end of the second production interval, on the twentyeighth day.The APB predicted in Annulus A was 65.0 MPa by the pseudo-steady-state model and 65.5 MPa by the hybrid Laplace model.In the other annuli, there is a slightly greater difference between the pressure increase predicted by the models, with 59. 34    To gain a better understanding of the differences between each model within shorter time frames, Table 3 presents the maximum APB difference observed in each annulus for the models.Additionally, it indicates the time at which this difference was detected.The table reveals that the speed of pressure buildup varies across annuli.Since the maximum pressure difference occurs at different times for each annulus, it is evident that the capacitance term slows down the pressure increase in the outer annuli.The hybrid Laplace model demonstrates an improvement in the results for small transients, which aligns with its main objective of reducing the overestimation of APB at the beginning of production.By incorporating the capacitance term in the heat equation for the annuli, pressurization should tend to zero within very short times, as there is insufficient time for heat to flow from the production tubing to the annulus and increase its temperature.This observation is confirmed for the outer annulus at 15 min of production, where the hybrid Laplace model predicts an APB of 2.80 MPa for Annulus B and 0.76 MPa for Annulus C.Although these values may still indicate significant pressurization, a comparison with the pseudo-steady-state model, which predicted 18.71 MPa and 14.65 MPa, respectively, clearly demonstrates that the new model does not suffer from the same high initial pressurization issue.

Conclusions
This paper presented a new formulation for modeling heat transfer at hydrocarbon producing wells based on the Laplace transformation of an equivalent thermal network with thermal capacity terms.This formulation was evaluated against real field data and compared in terms of physical coherence.From the analysis, the hybrid Laplace model was capable of accurately predicting the heat transfer in wellbores for long operations, with the capacity to reproduce the behavior of the traditional pseudo-steady-state formulation for long-time simulations.
When compared to field data in the vertical well, it was not possible to assert, with this study, if for a shorter timespan the hybrid Laplace model was more accurate.However, when analyzing the behavior of the APB in the annuli and the behavior of the temperature distribution in the wellbore, the hybrid Laplace predicted a behavior more aligned with the expected behavior dictated by the physics of the problem.
For the well in the EWT, the hybrid Laplace method was capable of better representing the thermal behavior of the produced fluid earlier than the traditional pseudo-steadystate formulation.Most of the deviations perceived in the data occurred in regions where there was lag or lack of data in the production flow rate.Therefore, improvements in this first analysis could involve a better representation of the oil flow rate.However, with the data provided, the hybrid model was better suited for simulating shorter operations.
For both cases, the APB increases in every annulus but at a slower rate when compared to the pseudo-steady-state.The maximum difference in APB is predicted at different times for each annulus, with the time of greatest difference increasing for the outer annuli.
A possible improvement for representing thermal behavior in even shorter operations would involve developing equations using the time-shifting property in Laplace transformation.By applying the transformation to delayed functions, it would be possible to incorporate phenomena such as varying oil flow rates into the equations and enhance the overall accuracy of the model, especially at the start of production where higher flow rates with larger variations are expected.This research is now under way and shall be presented in future publications.

Fig. 1
Fig. 1 Schematic diagram of a generic well with two annuli filled with packer fluids and one annulus cemented up to the wellhead.r 1 to r 7 are the radii of the interfaces between different materials/domains.The rectangular region (red contour) near the bottom of the well is used to illustrate the derivation of the governing equations

Fig.
Fig. Vertical well geometry

Fig. 5
Fig. 5 TPT data since start of production

Fig. 8
Fig. 8 Convergence test for Laplace model in the Vertical well (production time of 6 h)

Fig. 9 Fig. 10
Fig.9APB predicted with pseudo-steady-state and hybrid Laplace for vertical well

Fig. 11
Fig. 11 Temperature profile predicted by heat transfer models for 6 h of production

Fig. 12
Fig.12 Temperature profiles in the vertical well after 10 days of production

Fig. 13
Fig.13 Strain profile predicted by both models for vertical well at 6 h of production

Fig. 14 Fig. 15 Fig. 16
Fig.14 Strain profiles predicted by both models for vertical well at 10 days of production

Fig. 19
Fig. 19 APB predicted at each annuli for the EWT well by both proposed models

Table 1
Average temperature difference between the pseudo-steadystate and hybrid Laplace models for each annulus at different times

Table 2
Absolute average deviation of temperature and pressure measurements at the TPT across the production intervals MPa for the pseudo-steady-state model and 60.0 MPa for the hybrid Laplace model in Annulus B, and 56.13 MPa for the pseudo-steady-state model and 57.14 MPa for the hybrid Laplace model in Annulus C.

Table 3
Maximum difference between the APB predicted by the models for each annulus, and the time at which the maximum difference occurred