Study of the Thermobaric Model of Drilling Fluid in the Annular Space of Wells

1. Introduction

Deep and ultra-deep drilling, as well as the development of high-temperature and high-pressure fields, are key directions in the oil and gas industry. The efficiency and safety of drilling critically depend on the properties of the drilling fluid – a multifunctional dispersed system that performs a number of vital tasks: cuttings transport, maintaining hydrostatic pressure in the well, cooling the bit, and stabilizing the borehole wall [1,2,3,4].

During circulation, the drilling fluid is subjected to extreme thermobaric loads [4,5,6,7]. With increasing depth, temperature and pressure rise significantly, inevitably leading to changes in its physicochemical and structural-mechanical properties [8,9,10,11,12]. An increase in temperature can initiate thermal decomposition of components, coagulation, polymer degradation and, consequently, a sharp decrease in viscosity and rheological stability. At the same time, an increase in pressure usually causes compression of the dispersed phase and a reduction in interparticle distance, which, conversely, can lead to an increase in viscosity. These oppositely directed effects of $T$ and $P$ interact in a complex nonlinear manner. Accurate prediction of rheological parameters under HP/HT (high-pressure/high-temperature) conditions is a critical task for drilling control [12,13,14,15,16].

Field and laboratory evidence consistently shows that these effects cannot be neglected. Downhole pressure and temperature measurements in an offshore well drilled with a synthetic-based fluid demonstrated that the static and circulating density profiles deviate markedly from the surface-measured mud weight, because the fluid is heated and compressed simultaneously along the wellbore [5]. Systematic laboratory studies of base oils, brines and complete drilling fluids at extreme temperatures and pressures confirmed that the density of the liquid phase can change by several percent over the ranges typical of deep wells, with oil-based systems being noticeably more compressible and more thermally expansible than water-based ones [17,18]. Compositional models, in which the density of the whole fluid is assembled from the volumetric behaviour of its liquid components while the solids are treated as practically incompressible, have therefore become a common engineering approach [19]. More recently, automated laboratory set-ups have been developed to measure the pressure and temperature dependence of drilling-fluid density directly [20], and dedicated density–pressure–temperature formulations for the aqueous phase have been proposed for real-time HP/HT applications [21]. The reference equation of state for water [22] provides a rigorous thermodynamic basis for such formulations, but it is too cumbersome to be embedded directly in closed-form hydraulic calculations.

In parallel, a large body of work has addressed the prediction of downhole density with data-driven methods, including support vector machines, radial basis function and feed-forward neural networks, gradient-boosted trees and hybrid algorithms optimised by metaheuristics [23,24,25,26,27]. These models reproduce laboratory datasets with high accuracy, but their coefficients have no direct physical meaning, their reliability outside the range of the training data is uncertain, and they do not by themselves ensure consistency between density, pressure and flow regime along the wellbore. A physically transparent model with a small number of measurable parameters therefore remains of practical interest, both as a stand-alone design tool and as a physical constraint for data-driven approaches.

Heat and mass transfer processes, as well as hydraulic losses, are most intense in the annular space – the narrow annular gap between the drill string and the borehole wall (or casing) [28,29,30,31,32]. Here, the greatest increase in fluid temperature occurs, and hydraulic resistance (pressure drop) is formed, which directly affects such parameters as: cuttings removal efficiency, the magnitude of hydrodynamic pressures affecting borehole stability and preventing fluid influx, and circulation energy costs [33,34,35].

Incorrect accounting or inaccurate prediction of rheological and hydraulic parameters in the annulus under HP/HT conditions can lead to: reduced rate of penetration, difficulties in cuttings transport, increased accident rate (sticking, caving, lost circulation), and critical errors in cementing program design.

These consequences are most severe in wells with a narrow drilling window, where the admissible range of the equivalent circulating density (ECD) between the pore (or collapse) pressure gradient and the fracture pressure gradient may amount to only a few tens of kilograms per cubic metre. Under such conditions, an error in the downhole density or in the annular frictional pressure loss of the same order as the thermobaric correction itself can shift the bottomhole pressure outside the safe interval [34,36]. Real-time ECD estimation has therefore attracted considerable attention, ranging from hydraulic models calibrated against distributed downhole pressure data [35] to machine-learning estimators trained on surface drilling parameters [37].

Early works in the field of modelling rheological properties of drilling fluids focused on simple models: Bingham-Plastic and Power-Law. Currently, the generally accepted standard for most modern drilling fluids (especially oil-based ones) is the Herschel-Bulkley model [38]:

$$\tau ={\tau }_y+K{\dot{\gamma }}^n.$$

(1)

Historically, simple rheological models were widely used to describe drilling fluid flow. The Bingham model assumes a shear stress that must be overcome to initiate flow. While the Bingham model describes behaviour at high shear rates well, it is often inadequate at low shear rates. The Power-Law model, on the contrary, describes low shear rates better but does not have an explicit yield stress.

In modern practice, especially when working with polymer-based and oil-based muds, the Herschel–Bulkley (HB) model has become the gold standard, combining the advantages of both [2,38,39,40,41,42,43,44]. Using three parameters, this model provides the best approximation of experimental data over a wide range of shear rates. However, the true complexity lies not in the choice of the model, but in accurate determination how these three key parameters change under the extreme thermobaric conditions of deep wells.

The temperature and pressure dependence of the rheological parameters has been investigated experimentally for water-based, oil-based, formate and geothermal drilling fluids [16,45,46]. As a rule, the yield stress and the consistency index decrease with increasing temperature and increase with pressure, whereas the flow behaviour index changes only moderately; the magnitude of these trends, however, depends strongly on the fluid composition. The measured parameters are also sensitive to the viscometer geometry, and shear-rate corrections are required when Herschel–Bulkley fluids are tested in Couette viscometers [47]. Alternative multi-parameter rheological models, such as the vom Berg and Hahn–Eyring models, have been proposed to describe specific fluid systems more accurately [48,49], and in situ thermal conditions have been shown to affect surge and swab pressures significantly [50].

Several authors have proposed models that include both factors, often using multiplicative corrections. A typical approach consists of separate calculation of temperature and pressure factors [51,52,53,54,55,56]:

$$\mathsf{\mu }\left(\mathrm{P},\mathrm{T}\right)={\mathsf{\mu }}_{\mathrm{base}}\cdot \mathrm{f}\left(\mathrm{T}\right)\cdot \mathrm{g}\left(\mathrm{P}\right).$$

(2)

where $\mathrm{f}(\mathrm{T})$ is the thermal factor and $\mathrm{g}(\mathrm{P})$ is the baric factor.

Most empirical and semi-empirical models that account for the influence of $\mathrm{T}$ and $\mathrm{P}$ on rheological parameters are either highly specialised (applicable only to a specific fluid type or $\mathrm{T}/\mathrm{P}$ range) or use linear or additive dependencies, which do not reflect the complex mechanism of structural changes. Comprehensive mathematical models capable of adequately describing the nonlinear interdependence of rheological parameters on simultaneously acting $\mathrm{T}$ and $\mathrm{P}$ are insufficiently developed, especially at high shear rates characteristic of the annulus. Hydraulic calculations in the annulus require integration of the equations of motion in the annular gap, which is complicated by the following factors [57,58,59,60,61]: non-Newtonian flow behaviour (use of the Herschel–Bulkley model), annular eccentricity (non--concentric arrangement of strings), non-isothermal flow (significant temperature gradient along the annulus), and mixed flow regime (transition from laminar to turbulent flow).

Existing hydraulic models often neglect a detailed account of changes in rheological parameters along the annulus due to temperature and pressure effects, which leads to significant errors in pressure loss calculations and, consequently, in the estimation of bottomhole hydrodynamic pressures.

A further motivation for compact thermobaric models comes from the rapid development of automated drilling-fluid monitoring. Inline and skid-based systems that measure density, rheology and other properties at the rig site at intervals of several minutes are now available and have been tested in field conditions [62,63], including HP/HT gas wells drilled with managed pressure drilling (MPD) systems [64,65]. The authors’ previous work has addressed measurement devices, improved viscometric procedures and software–hardware complexes for continuous monitoring of drilling-fluid properties [31,32,66,67,68]. However, all such systems measure the fluid at surface conditions. To use these data for downhole pressure management, the measured values must be converted to the actual temperature and pressure at respective depth, which requires a model that is both physically sound and computationally inexpensive.

The annular space is a critical modelling object. It is the zone, where temperature gradients are most pronounced and where bottomhole pressure has the greatest impact on well stability. The drilling fluid is the most important component of the drilling system, responsible for cuttings transport, borehole wall stabilisation and, most importantly, for bottomhole pressure control. When drilling deep wells or wells with high pressure and high temperature, traditional surface measurements of rheology and density lose their accuracy. The plastic viscosity and density can change significantly under the pressure (compaction) and temperature (thermal expansion). Inaccurate measurement of rheological parameters in the annulus leads to erroneous pressure loss calculations, thereby narrowing the drilling window. Thus, the aim of this work is to develop a mathematical model for accurate prediction of rheological and hydraulic parameters in the annulus as functions of depth $\mathrm{H}$.

The novelty of the present study covers the following: (i) the linearity of the thermal expansion coefficient with respect to temperature and of the isothermal bulk modulus with respect to pressure is verified statistically for the range of conditions typical of wells deeper than 1000 m; (ii) on this basis, the differential equation of volume change is integrated in closed form along the annulus, with the local hydraulic resistance included separately for laminar and turbulent flow; and (iii) the resulting density equations are combined with the drilling-window constraints into systems that can be solved directly for the admissible circulation velocity and fluid parameters.

2. Materials and Methods

The coefficient α of thermal expansion (with respect to temperature change), ${\mathrm{K}}^{-1}$, and the isothermal compressibility coefficient β (with respect to pressure change), ${\mathrm{MPa}}^{-1}$, characterize the change in fluid volume. Thus, the change in fluid volume $\mathrm{V}$ depends on temperature $\mathrm{T}$ and pressure $\mathrm{P}$:

V=V(P,T).

(3)

From mathematical analysis, the total differential of the function of two variables describes its change:

$$\mathrm{d}\mathrm{V}={\displaystyle \frac{\partial \mathrm{V}(\mathrm{P},\mathrm{T})}{\partial \mathrm{T}}}\mathrm{d}\mathrm{T}+{\displaystyle \frac{\partial \mathrm{V}(\mathrm{P},\mathrm{T})}{\partial \mathrm{P}}}\mathrm{d}\mathrm{P}.$$

(4)

Dividing by $\mathrm{V}$ gives:

$${\displaystyle \frac{\mathrm{d}\mathrm{V}}{\mathrm{V}}}={\displaystyle \frac{1}{\mathrm{V}}}\cdot {\displaystyle \frac{\partial \mathrm{V}(\mathrm{P},\mathrm{T})}{\partial \mathrm{T}}}\mathrm{d}\mathrm{T}+{\displaystyle \frac{1}{\mathrm{V}}}\cdot {\displaystyle \frac{\partial \mathrm{V}(\mathrm{P},\mathrm{T})}{\partial \mathrm{P}}}\mathrm{d}\mathrm{P}.$$

(5)

Here, ${\displaystyle \frac{1}{\mathrm{V}}}\cdot {\displaystyle \frac{\partial \mathrm{V}(\mathrm{P},\mathrm{T})}{\partial \mathrm{T}}}=\mathsf{\alpha }$ is the coefficient of thermal expansion, and $-{\displaystyle \frac{1}{\mathrm{V}}}\cdot {\displaystyle \frac{\partial \mathrm{V}(\mathrm{P},\mathrm{T})}{\partial \mathrm{P}}}=\mathsf{\beta }$ is the isothermal compressibility coefficient. Hence, it can be written as follows:

$${\displaystyle \frac{\mathrm{d}\mathrm{V}}{\mathrm{V}}}=\mathsf{\alpha }\hspace{0.17em}\mathrm{d}\mathrm{T}-\mathsf{\beta }\hspace{0.17em}\mathrm{d}\mathrm{P}.$$

(6)

In wells deeper than 1000 m, the temperature exceeds 20 °C and pressure exceeds 10 MPa. Therefore, the functions $\mathsf{\alpha }=\mathsf{\alpha }(\mathrm{T})$ and $\mathsf{\beta }=\mathsf{\beta }(\mathrm{P})$ must be determined. First, we evaluate these dependencies for linearity. If linear, we derive the explicit relations. We consider $\mathsf{\alpha }=\mathsf{\alpha }(\mathrm{T})$ at constant $\mathrm{P}$ and $\mathsf{\beta }=\mathsf{\beta }(\mathrm{P})$ at constant $\mathrm{T}$. The results are combined with Eq. (6). The linearity assessment was performed using an experimental rig shown in Figure 1, filled with wellbore water, and the conclusions are considered valid for water-based drilling fluids.

2.1. Experimental Setup and Measurement Procedure

The laboratory rig shown in Figure 1 reproduces the thermobaric conditions of the annular space. The wellbore simulator body (6), closed by the top cover (17), contains a drill-string model (7, 20, 21) suspended on a rod (16) and surrounded by the test fluid (8), so that an annular gap is formed between the string and the simulator wall. The temperature is controlled by the circulating temperature-control liquid (3), which is heated by the electric heater (5) and circulated by the pump (12) driven by the motor (13) through the control valve (14); the thermal insulation pad (2) reduces heat losses to the rig base. The pressure is generated by the hydraulic press (9) connected to the simulator through the high-pressure tube (9a) and the press coupling (15). The tube (11) connecting the tank (4) with the simulator makes it possible to reproduce water influx. Temperature and pressure are measured by the thermometer (19) and the pressure gauge (22); all signals are collected by the microcontroller (25) and transmitted via the communication cable (24) to the process-control computer (10). The rig allows tests at temperatures from 1 to 110 °C and at pressures up to 50 MPa; the thermal expansion series were run at 0.1, 10, 20 and 50 MPa and the compressibility series at 0.5–8 MPa, as indicated in Table 1 and Table 3.

The thermal expansion coefficient was determined at constant pressure by raising the temperature from 0 to 100 °C in steps of 10 °C. The series were repeated at pressures of 0.1, 10, 20 and 50 MPa, and results are collected in Table 1 and discussed in the next section. The isothermal compressibility was determined at constant temperature by raising the pressure in the steps 0.5, 1, 2, 4 and 8 MPa. These series were repeated at temperatures of 0, 5, 10, 15 and 25 °C. After each step, the readings were held for 15 min until the temperature and the pressure had stabilised, and the corresponding change in the fluid volume was recorded. Each series was repeated three times, and the arithmetic mean values were used in the correlation analysis. The temperature was measured with a laboratory thermometer with a scale interval of 1 °C and a maximal error of ±0.5 °C. For the pressure, a gauge of accuracy class 0.6 was used with an upper measuring limit of 60 MPa, i.e. with a limit of error of ±0.36 MPa.

The density of the test fluid at surface conditions, ρ₀, was determined with an ABR-1 hydrometer. The permissible additional error of this instrument does not exceed 2 kg/m³ for every 10 °C change of the test-solution temperature from 20±2 °C, and 10 kg/m³ for every 10 °C change caused by climatic environmental factors. The hydrometer was calibrated with distilled water at 20±5 °C, for which the reading must equal 1000 kg/m³. These uncertainties propagate directly into ρ₀ and, therefore, into the density profiles calculated with the model.

The linearity analysis was performed with wellbore water. This choice is justified by the fact that, in water-based drilling fluids, the liquid phase determines the volumetric response to temperature and pressure. The solid phase consisting of clay and weighting material is practically incompressible and has a thermal expansion coefficient about an order of magnitude lower than that of water. The same principle underlies compositional density models [17,19].

2.2. Linear Dependence of the Thermal Expansion Coefficient α=α(T) at Constant P_hydr for Wellbore Water

Experimental values of the thermal expansion coefficient α×${10}^{-4}\hspace{0.17em}{\mathrm{K}}^{-1}$ for wellbore water at different temperatures and pressures are given in Table 1. The measured values of Table 1 agree with the reference equation of state for water [22] within about 1% over most of the range investigated. For oil-based and synthetic fluids, whose base liquids are considerably more compressible, the parameters must be determined separately (see Section 4.4).

To evaluate linearity of $\mathsf{\alpha }=\mathsf{\alpha }(\mathrm{T})$, the correlation coefficient is calculated as:

$${\mathrm{r}}_c={\displaystyle \frac{\mathrm{n}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{x}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}-\left({\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{x}}_{\mathrm{i}}\right)\left({\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{y}}_{\mathrm{i}}\right)}{\sqrt{\mathrm{n}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{x}}_{\mathrm{i}}^2-{\left({\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{x}}_{\mathrm{i}}\right)}^2}\hspace{0.25em} \sqrt{\mathrm{n}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{y}}_{\mathrm{i}}^2-{\left({\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{y}}_{\mathrm{i}}\right)}^2}}}.$$

(7)

The correlation coefficient lies in the range $\left[-1,1\right]$, and its absolute value can be interpreted as follows:

$0.99\le ⌊{\mathrm{r}}_c⌋\le 1$ – excellent correlation,

$0.95\le ⌊{\mathrm{r}}_c⌋<0.99$ – good correlation,

$0.9<⌊{\mathrm{r}}_c⌋<0.95$ – satisfactory correlation,

$\mid {\mathrm{r}}_c\mid \le 0.9$ – no linear relationship.

For high accuracy, the least-squares method is recommended only when $0.985\le ⌊{\mathrm{r}}_c⌋\le 1$. Using Table 1 and Eq. (7), the correlation coefficients and linearity assessments are as follows:

Table 2. Correlation coefficients for $\mathsf{\alpha }=\mathsf{\alpha }(\mathrm{T})$.

Pressure $\mathbf{P}$, MPa	rc	Assessment
0.1	0.9803	NOT LINEAR
10	0.9899	LINEAR
20	0.9906	LINEAR
50	0.994	LINEAR

Table 2. Correlation coefficients for $\mathsf{\alpha }=\mathsf{\alpha }(\mathrm{T})$.

Pressure $\mathbf{P}$, MPa	rc	Assessment
0.1	0.9803	NOT LINEAR
10	0.9899	LINEAR
20	0.9906	LINEAR
50	0.994	LINEAR

Higher pressure reduces the error of the linear relationship between temperature and thermal expansion coefficient. In wells deeper than 1000 m, pressure exceeds 10 MPa and temperature exceeds 20 °C, so $\mathsf{\alpha }=\mathsf{\alpha }(\mathrm{T})$ is linear for wellbore water with an error less than 1 %.

For a linear dependence, the derivative is constant. Since the above evaluation shows that $\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})$ is constant for wellbore water in deep wells:

$$\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})={\displaystyle \frac{\mathrm{d}\mathsf{\alpha }(\mathrm{T})}{\mathrm{d}\mathrm{T}}}=\mathrm{const}.$$

(8)

Here $\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})$ (in ${\mathrm{K}}^{-2}$) is the gradient of the thermal expansion coefficient of the drilling fluid at a given pressure, determined in the laboratory from:

$${\displaystyle \frac{\mathsf{\Delta }\mathsf{\alpha }(\mathrm{T})}{\mathsf{\Delta }\mathrm{T}}}={\displaystyle \frac{\mathrm{d}\mathsf{\alpha }(\mathrm{T})}{\mathrm{d}\mathrm{T}}}=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T}).$$

(9)

$\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})$ is constant for a given drilling fluid. Integrating (8) gives:

$$\mathsf{\alpha }(\mathrm{T})=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})\cdot \mathrm{T}+\mathsf{\alpha }({\mathrm{T}}_0,{\mathrm{P}}_0).$$

(10)

Equations (8–10) are valid only for temperatures above 20 °C and pressures above 10 MPa, i.e., for well depths greater than 1000 m.

2.3. Linear Dependence of the Isothermal Compressibility Coefficient

Empirical values of the isothermal compressibility coefficient β, ${\mathrm{MPa}}^{-1}$, for wellbore water are given in Table 3. These measurements were performed at pressures from 0.5 MPa to 8 MPa and temperatures from 0 °C to 25 °C, whereas the thermal expansion data of Table 1 cover pressures up to 50 MPa and temperatures up to 100 °C. The extrapolation of the resulting linear law to the pressures encountered in deep wells.

Using Eq. (7), the correlation coefficients r_c for $\mathsf{\beta }(\mathrm{P})$ are calculated and collected in Table 4. However, better linear correlation is found for the isothermal bulk modulus ${\mathrm{E}}_{\mathrm{up}}=1/\mathsf{\beta }$ as a function of pressure, as shown in Table 5.

Within the investigated range, a higher temperature reduces the error of the linear relationship between pressure and the isothermal bulk modulus: at temperatures of 5 °C and above, the linearity criterion (7) is satisfied, whereas at 0 °C it is not. Since the temperature of a well deeper than 1000 m exceeds 20 °C, the linear approximation is adopted in the form:

$$\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}({\mathrm{E}}_{\mathrm{up}})={\displaystyle \frac{\mathrm{d}({\mathrm{E}}_{\mathrm{up}})}{\mathrm{d}\mathrm{P}}}=\mathrm{const},$$

(11)

where ${\mathrm{E}}_{\mathrm{up}}=1/\mathsf{\beta }$, MPa. Then:

$$\mathrm{d}({\mathrm{E}}_{\mathrm{up}})=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}({\mathrm{E}}_{\mathrm{up}})\hspace{0.17em}\mathrm{d}\mathrm{P}.$$

(12)

Integrating:

$${\mathrm{E}}_{\mathrm{up}}=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}({\mathrm{E}}_{\mathrm{up}})\cdot \mathrm{P}+{\mathrm{E}}_{\mathrm{up}}({\mathrm{P}}_0,{\mathrm{T}}_0).$$

(13)

Hence:

$$\mathsf{\beta }(\mathrm{P})={\displaystyle \frac{1}{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}({\mathrm{E}}_{\mathrm{up}})\cdot \mathrm{P}+{\mathrm{E}}_{\mathrm{up}}({\mathrm{P}}_0,{\mathrm{T}}_0)}}.$$

(14)

Equations above are valid for temperatures T > 20 °C and pressure $\mathrm{P}>10\hspace{0.17em}\mathrm{MPa}$, corresponding to wells deeper than 1000 m.

2.4. Pressure Balance During Circulation

The bottomhole pressure during circulation is given by:

$${\mathrm{P}}_{\mathrm{b}}={\mathrm{P}}_{\mathrm{h}}+\mathsf{\Delta }{\mathrm{P}}_{\mathrm{f}}+{\mathrm{P}}_0,$$

(15)

where ${\mathrm{P}}_{\mathrm{h}}$ is the hydrostatic pressure (MPa), $\mathsf{\Delta }{\mathrm{P}}_{\mathrm{f}}$ the frictional pressure loss in the annulus (MPa), and ${\mathrm{P}}_0$ the surface pressure (MPa). Because the fluid moves upward in the annulus, the frictional losses increase the bottomhole pressure, so that the equivalent circulating density always exceeds the static one; this sign convention is used consistently in Equations (23), (29)–(36). The hydrostatic pressure is:

$${\mathrm{P}}_{\mathrm{h}}=\mathsf{\rho }\mathrm{g}\mathrm{H},$$

(16)

with $\mathsf{\rho }$ the drilling fluid density, kg/m³; $\mathrm{g}$ gravity acceleration, m/s²; and $\mathrm{H}$ the well depth, m. The frictional loss is described by the Darcy-Weisbach equation:

$$\mathsf{\Delta }{\mathrm{P}}_{\mathrm{f}}={\displaystyle \frac{\mathsf{\lambda }\mathsf{\rho }{\mathrm{v}}^2\mathrm{H}}{2{\mathrm{d}}_{\mathrm{h}}}}.$$

(17)

Here $\mathsf{\lambda }$ is the friction factor; $\mathrm{v}$ the mean flow velocity, m/s; and ${\mathrm{d}}_{\mathrm{h}}={\mathrm{d}}_{\mathrm{w}}-{\mathrm{d}}_{\mathrm{o}}$ the hydraulic diameter of the annulus, m; with ${\mathrm{d}}_{\mathrm{w}}$ the wellbore diameter and ${\mathrm{d}}_{\mathrm{o}}$ the outer drill pipe diameter, respectively. The Reynolds number in the annulus is:

$$\mathrm{R}\mathrm{e}={\displaystyle \frac{\mathsf{\rho }\mathrm{v}{\mathrm{d}}_{\mathrm{h}}}{\mathsf{\mu }}}={\displaystyle \frac{\mathrm{v}{\mathrm{d}}_{\mathrm{h}}}{\mathsf{\nu }}},$$

(18)

where $\mathsf{\mu }$ is the dynamic viscosity, Pa·s; and $\mathsf{\nu }=\mathsf{\mu }/\mathsf{\rho }$ is the kinematic viscosity, m²/s.

For laminar flow ($\mathrm{R}\mathrm{e}\le 2300$), $\mathsf{\lambda }=64/\mathrm{R}\mathrm{e}$, and Eq. (17) becomes:

$$\mathsf{\Delta }{\mathrm{P}}_{\mathrm{f}}={\displaystyle \frac{32\mathsf{\mu }\mathrm{v}\mathrm{H}}{{\mathrm{d}}_{\mathrm{h}}^2}}.$$

(19)

For turbulent flow in smooth pipes ($\mathrm{R}\mathrm{e}>4000$), the Blasius approximation gives $\mathsf{\lambda }=0.3164/\mathrm{R}{\mathrm{e}}^{1/4}$. Substituting into Eq. (17):

$$\mathsf{\Delta }{\mathrm{P}}_{\mathrm{f}}={\displaystyle \frac{0.1582\hspace{0.17em}{\mathsf{\rho }}^{3/4}\hspace{0.17em}{\mathsf{\mu }}^{1/4}\hspace{0.17em}{\mathrm{v}}^{7/4}\mathrm{H}}{{\mathrm{d}}_{\mathrm{h}}^{5/4}}}.$$

(20)

2.5. Temperature and Pressure Profiles

The geothermal gradient is assumed constant:

$${\displaystyle \frac{\mathrm{d}\mathrm{T}}{\mathrm{d}\mathrm{H}}}=\mathrm{grad}\mathrm{T}=\mathrm{m}=\mathrm{const}$$

(21)

Thus, the temperature at depth $\mathrm{H}$ is:

$$\mathrm{T}=\mathrm{m}\mathrm{H}+{\mathrm{T}}_0,$$

(22)

where ${\mathrm{T}}_0$ is the surface temperature, K. From the pressure balance (15)–(17), the pressure at depth $\mathrm{H}$ can be expressed as:

$$\mathrm{P}=\mathsf{\xi }(\mathrm{H})\mathrm{H}+{\mathrm{P}}_0,\mathsf{\xi }(\mathrm{H})=\mathsf{\rho }\mathrm{g}+{\displaystyle \frac{\mathsf{\lambda }\mathsf{\rho }{\mathrm{v}}^2}{2{\mathrm{d}}_{\mathrm{h}}}}.$$

(23)

After differentiating, it becomes as follows:

$$\mathrm{d}\mathrm{P}=(\mathrm{H}{\mathsf{\xi }}^{\prime }(\mathrm{H})+\mathsf{\xi }(\mathrm{H}))\mathrm{d}\mathrm{H},\mathrm{d}\mathrm{T}=\mathrm{m}\hspace{0.17em}\mathrm{d}\mathrm{H}.$$

(24)

2.6. Density Model for Circulating Conditions

Substituting (24) into the differential volume change Eq. (6), it appears as follows:

$${\displaystyle \frac{\mathrm{d}\mathrm{V}}{\mathrm{V}}}=\mathsf{\alpha }\hspace{0.17em}\mathrm{d}\mathrm{T}-\mathsf{\beta }\hspace{0.17em}\mathrm{d}\mathrm{P}=\mathsf{\alpha }\hspace{0.17em}\mathrm{m}\hspace{0.17em}\mathrm{d}\mathrm{H}-\mathsf{\beta }(\mathrm{H}{\mathsf{\xi }}^{\prime }(\mathrm{H})+\mathsf{\xi }(\mathrm{H}))\mathrm{d}\mathrm{H}.$$

(25)

Using the linear relations (10) for α and (13) for Eup $\mathsf{\alpha }(\mathrm{T})\mathsf{\beta }(\mathrm{P})$, and after algebraic manipulation (see Appendix A), the density at depth $\mathrm{H}$ is obtained as:

$$\mathsf{\rho }(\mathrm{H})={\displaystyle \frac{{\mathsf{\rho }}_0\hspace{0.17em}\sqrt[\mathrm{q}]{\frac{\mathrm{q}\hspace{0.17em}\mathsf{\xi }(\mathrm{H})\mathrm{H}}{\mathrm{q}{\mathrm{P}}_0+\mathrm{r}}+1}}{\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{{\mathrm{m}}^2\mathrm{n}{\mathrm{H}}^2}{2}+\mathrm{m}\mathrm{n}{\mathrm{T}}_0\mathrm{H}+\mathrm{m}\mathrm{p}\mathrm{H}\right)}},$$

(26)

where the constants are defined as:

$$\mathrm{m}=\mathrm{grad}\mathrm{T},\mathrm{n}=\mathrm{grad}\mathsf{\alpha }\left(\mathrm{T}\right),\mathrm{p}=\mathsf{\alpha }\left({\mathrm{T}}_0,{\mathrm{P}}_0\right),\mathrm{q}=\mathrm{grad}\left({\mathrm{E}}_{\mathrm{up}}\right),\mathrm{r}={\mathrm{E}}_{\mathrm{up}}\left({\mathrm{P}}_0,{\mathrm{T}}_0\right),$$

(27)

and ${\mathsf{\rho }}_0$ is the fluid density at surface conditions (${\mathrm{T}}_0,{\mathrm{P}}_0$).

2.6.1. Laminar Flow Regime

Insert the laminar expression for ξ derived from (19):

$$\mathsf{\xi }\left(\mathrm{H}\right)=\mathsf{\rho }\mathrm{g}+{\displaystyle \frac{32\mathsf{\mu }\mathrm{v}}{{\mathrm{d}}_{\mathrm{h}}^2}},$$

(28)

with $\mathsf{\mu }=\mathrm{A}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{E}/(\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})))$ from the Frenkel-Andrade equation, where A is the pre-exponential factor, Pa·s; E is the activation energy of viscous flow, and k is the Boltzmann constant (or the universal gas constant when E is expressed per mole). Substituting into (26) and raising to the power q yields the algebraic equation:

$$\left\{\begin{array}{c}\mathrm{b}{\mathsf{\rho }}^{\mathrm{q}}-\mathrm{\rho }+\mathrm{a}=0\\ \mathrm{\rho }\mathrm{v}\le {\displaystyle \frac{2300\hspace{0.17em}\mathrm{A}\mathrm{e}\mathrm{x}\mathrm{p}(\frac{\mathrm{E}}{\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})})}{{\mathrm{d}}_{\mathrm{h}}}}\end{array}\right..$$

(29)

with

$$\mathrm{a}=-{\displaystyle \frac{32\mathrm{A}\mathrm{v}}{{\mathrm{d}}_{\mathrm{h}}^2\mathrm{g}}}\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{\mathrm{E}}{\mathrm{k}\left({\mathrm{T}}_0+\mathrm{m}\mathrm{H}\right)}}\right)-{\displaystyle \frac{{\mathrm{P}}_0}{\mathrm{g}\mathrm{H}}}-{\displaystyle \frac{\mathrm{r}}{\mathrm{q}\mathrm{g}\mathrm{H}}},\mathrm{b}={\displaystyle \frac{\mathrm{q}{\mathrm{P}}_0+\mathrm{r}}{{\mathsf{\rho }}_0^{\mathrm{q}}\mathrm{q}\mathrm{H}\mathrm{g}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\mathrm{m}}^2\mathrm{n}\mathrm{q}{\mathrm{H}}^2}{2}}-\mathrm{m}\mathrm{n}{\mathrm{T}}_0\mathrm{q}\mathrm{H}-\mathrm{m}\mathrm{p}\mathrm{q}\mathrm{H}\right).$$

(30)

The inequality ensures laminar flow at $\mathrm{R}\mathrm{e}\le 2300$.

2.6.2. Turbulent Flow Regime

Using the turbulent expression (20) and $\mathsf{\mu }=\mathrm{A}\mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{E}/(\mathrm{k}\mathrm{T}))$, the hydraulic resistance term becomes:

$$\mathsf{\xi }(\mathrm{H})=\mathsf{\rho }\mathrm{g}+{\displaystyle \frac{0.1582\hspace{0.17em}{\mathsf{\rho }}^{3/4}(\mathrm{A}{\mathrm{e}}^{\mathrm{E}/(\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H}))}{)}^{1/4}{\mathrm{v}}^{7/4}}{{\mathrm{d}}_{\mathrm{h}}^{5/4}}}.$$

(31)

Inserting this into (26) and raising to the power q leads to following:

$$\left\{\begin{array}{c}\mathrm{c}{\mathsf{\rho }}^{\mathrm{q}}-\mathrm{\rho }-\mathrm{a}\sqrt[4]{{\mathsf{\rho }}^3}-\mathrm{b}=0\\ \mathrm{\rho }\mathrm{v}>{\displaystyle \frac{4000\hspace{0.17em}\mathrm{A}\mathrm{e}\mathrm{x}\mathrm{p}(\frac{\mathrm{E}}{\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})})}{{\mathrm{d}}_{\mathrm{h}}}}\end{array}\right..$$

(32)

where

$$\mathrm{a}={\displaystyle \frac{0.1582\sqrt[4]{\mathrm{A}}}{\mathrm{g}\sqrt[4]{{\mathrm{v}}^7{\mathrm{d}}_{\mathrm{h}}^5}}}\mathrm{e}\mathrm{x}\mathrm{p}({\displaystyle \frac{\mathrm{E}}{4\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})}}),\mathrm{b}={\displaystyle \frac{{\mathrm{P}}_0}{\mathrm{g}\mathrm{H}}}+{\displaystyle \frac{\mathrm{r}}{\mathrm{q}\mathrm{g}\mathrm{H}}},\mathrm{c}={\displaystyle \frac{\mathrm{q}{\mathrm{P}}_0+\mathrm{r}}{{\mathsf{\rho }}_0^{\mathrm{q}}\mathrm{g}\mathrm{q}\mathrm{H}}}\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \frac{{\mathrm{m}}^2\mathrm{n}\mathrm{q}{\mathrm{H}}^2}{2}}-\mathrm{m}\mathrm{n}{\mathrm{T}}_0\mathrm{q}\mathrm{H}-\mathrm{m}\mathrm{p}\mathrm{q}\mathrm{H}).$$

(33)

The inequality enforces turbulent flow at $\mathrm{R}\mathrm{e}>4000$.

2.7. Drilling Window and Safe Operating Conditions

To avoid formation fracture (lost circulation) or collapse (influx), the bottomhole pressure must satisfy the following:

$${\mathrm{P}}_{\mathrm{form}}<{\mathrm{P}}_{\mathrm{b}}<{\mathrm{P}}_{\mathrm{frac}},$$

(34)

where ${\mathrm{P}}_{\mathrm{form}}={\mathrm{K}}_{\mathrm{an}}{\mathsf{\rho }}_{\mathrm{w}}\mathrm{g}\mathrm{H}$, with ${\mathrm{K}}_{\mathrm{an}}$ the abnormality coefficient and ${\mathsf{\rho }}_{\mathrm{w}}$ water density, while ${\mathrm{P}}_{\mathrm{frac}}=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{\mathrm{P}}_{\mathrm{frac}}\cdot \mathrm{H}$. Combining this with the appropriate flow regime equation gives the final coupled systems.

For laminar flow:

$$\left\{\begin{array}{c}{\mathrm{K}}_{\mathrm{an}}{\mathsf{\rho }}_{\mathrm{w}}\mathrm{g}<\mathrm{\rho }\mathrm{g}+{\displaystyle \frac{32\mathrm{A}\mathrm{v}}{{\mathrm{d}}_{\mathrm{h}}^2}}\mathrm{e}\mathrm{x}\mathrm{p}({\displaystyle \frac{\mathrm{E}}{\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})}})+{\displaystyle \frac{{\mathrm{P}}_0}{\mathrm{H}}}<\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{\mathrm{P}}_{\mathrm{frac}}\\ \mathrm{b}{\mathsf{\rho }}^{\mathrm{q}}-\mathrm{\rho }+\mathrm{a}=0\\ \mathrm{\rho }\mathrm{v}\le {\displaystyle \frac{2300\mathrm{A}}{{\mathrm{d}}_{\mathrm{h}}}}\mathrm{e}\mathrm{x}\mathrm{p}({\displaystyle \frac{\mathrm{E}}{\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})}})\end{array}\right..$$

(35)

For turbulent flow:

$$\left\{\begin{array}{c}{\mathrm{K}}_{\mathrm{an}}{\mathsf{\rho }}_{\mathrm{w}}\mathrm{g}<\mathrm{\rho }\mathrm{g}+{\displaystyle \frac{0.1582\hspace{0.17em}{\mathsf{\rho }}^{3/4}(\mathrm{A}{\mathrm{e}}^{\mathrm{E}/(\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H}))}{)}^{1/4}{\mathrm{v}}^{7/4}}{{\mathrm{d}}_{\mathrm{h}}^{5/4}}}+{\displaystyle \frac{{\mathrm{P}}_0}{\mathrm{H}}}<\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{\mathrm{P}}_{\mathrm{frac}}\\ \mathrm{c}{\mathsf{\rho }}^{\mathrm{q}}-\mathrm{\rho }-\mathrm{a}\sqrt[4]{{\mathsf{\rho }}^3}-\mathrm{b}=0\\ \mathrm{\rho }\mathrm{v}>{\displaystyle \frac{4000\mathrm{A}}{{\mathrm{d}}_{\mathrm{h}}}}\mathrm{e}\mathrm{x}\mathrm{p}({\displaystyle \frac{\mathrm{E}}{\mathrm{k}({\mathrm{T}}_0+\mathrm{m}\mathrm{H})}})\end{array}\right..$$

(36)

The constants in (29) and (32) are those defined for the respective regimes above. These systems allow one to select the flow velocity $\mathrm{a},\mathrm{b},\mathrm{c}\mathrm{v}$ and fluid parameters that maintain a safe downhole pressure while keeping the desired flow regime in the annulus.

All models are valid for T > 20 °C and $\mathrm{P}>10\hspace{0.17em}\mathrm{MPa}$, i.e., for well depths greater than 1000 m, with an overall error not exceeding ±1.5 % when the laboratory data are accurate.

2.8. Determination of Model Parameters

The model contains three groups of input parameters shown in Table 6. Geometric and operational parameters, i.e. the well depth H, the hydraulic diameter d_h, the mean annular velocity v and the surface pressure P₀, are known from the well design and the drilling programme. The thermal parameter m is the temperature gradient along the annulus, which can be estimated from offset-well data or from downhole temperature measurements [4]. The fluid parameters are determined in the laboratory. The coefficients gradα(T) and α(T₀, P₀) are obtained by least-squares fitting of the measured α–T data at the pressure of interest, and grad(E_up) and E_up(P₀, T₀) by fitting of the E_up–P data at the temperature of interest; in both cases, the linear approximation is accepted only if |r_c| ≥ 0.985, see Equation (7). The constants A and E of the Frenkel–Andrade equation are found from viscosity measurements at no fewer than three temperatures by linear regression of ln μ against 1/T, whose slope equals E/k and whose intercept equals ln A. Finally, ρ₀ is the density measured at surface conditions.

2.9. Numerical Solution Procedure

For a given depth, Equations (29) and (32) are nonlinear algebraic equations with respect to ρ, which do not admit an explicit solution for an arbitrary exponent q. They are solved numerically by the Newton–Raphson method, with the surface density ρ₀ used as the initial approximation; because the density changes by only a few percent along the wellbore, convergence to a relative tolerance of 10⁻⁶ is usually achieved within a few iterations. The calculation proceeds as follows: (1) the temperature at depth H is calculated from Equation (22) and the viscosity from the Frenkel–Andrade equation; (2) a flow regime is assumed, and the corresponding equation, (29) or (32), is solved for ρ; (3) the Reynolds-number inequality of the selected system is checked, and if it is violated, the calculation is repeated for the other regime; (4) the bottomhole pressure P_b and the equivalent circulating density ECD = P_b/(gH) are evaluated and compared with the drilling-window constraints of Systems (35) and (36); and (5) if the constraints are violated, the circulation velocity v or the surface pressure P₀ is adjusted and the procedure is repeated.

In the transitional interval 2300 < Re < 4000, neither Equation (19) nor Equation (20) is strictly applicable. In engineering calculations, it is advisable to evaluate both and to adopt the larger frictional pressure loss as a conservative estimate. Repeating the procedure over a grid of depths yields the profiles of density, viscosity and pressure along the annulus. Because only one scalar equation is solved at each depth, the computational cost is negligible, which makes the procedure suitable for real-time applications.

3. Results

3.1. Linearity of the Thermal Expansion Coefficient and the Bulk Modulus

The experimental data for wellbore water collected in Table 1, Table 2, Table 3, Table 4 and Table 5 show that the coefficient of thermal expansion $\mathsf{\alpha }$ is a linear function of temperature when the pressure exceeds 10 MPa. For pressures of 10, 20 and 50 MPa, the correlation coefficient ${\mathrm{r}}_c$ lies between 0.9899 and 0.994, well within the linearity criterion (${\mathrm{r}}_c\ge 0.985$). At atmospheric pressure 0.1 MPa the linearity is poorer, and ${\mathrm{r}}_c=0.9803$, confirming that the linear approximation becomes excellent only under elevated pressures typical of deep wells. The temperature gradient of the expansion coefficient, $\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})$, is therefore a constant for a given drilling fluid under HP/HT conditions, allowing the simple relation $\mathsf{\alpha }(\mathrm{T})=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\mathsf{\alpha }(\mathrm{T})\cdot \mathrm{T}+\mathsf{\alpha }({\mathrm{T}}_0,{\mathrm{P}}_0)$ with an error below 1 %.

For the isothermal bulk modulus ${\mathrm{E}}_{\mathrm{up}}=1/\mathsf{\beta }$, the correlation with pressure is linear at temperatures of 5 °C and above (Table 5). At the lowest temperature (0 °C) the linearity is marginal, but under realistic downhole conditions (T > 20 °C) the linear approximation is justified with an error ≤ 1.5 %. It should be emphasized, however, that the compressibility data of Table 3, Table 4 and Table 5 were obtained at pressures up to 8 MPa; the use of the linear law E_up(P) at the pressures of 10–50 MPa typical of deep wells is therefore an extrapolation, whose experimental verification over the full pressure range of the rig is a subject of further work. This leads to $E_{\mathrm{up}}(P)=\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}(E_{\mathrm{up}})\cdot P+E_{\mathrm{up}}(P_0,T_0)$ and consequently to the pressure-dependent compressibility. These linearizations are a key simplification that makes the subsequent integration tractable while maintaining engineering accuracy.

From a physical standpoint, the poorer linearity of α(T) at atmospheric pressure is explained by the anomalous behaviour of water at low temperatures: near the temperature of maximum density (about 4 °C at 0.1 MPa), the thermal expansion coefficient changes sign and varies strongly with temperature. Increasing pressure shifts the temperature of maximum density to lower values and weakens this anomaly [22], so that in the temperature range typical of deep wells, α(T) becomes nearly linear. Similarly, the approximately linear increase in the bulk modulus with pressure is consistent with the Murnaghan- and Tait-type equations of state commonly used for liquids, in which the tangent bulk modulus is approximately a linear function of pressure. These physical arguments support the extension of the statistical conclusions obtained for wellbore water to the aqueous phase of water-based drilling fluids.

3.2. Density Prediction Under Circulation

The central result of this work is the closed-form expression for the drilling fluid density at any depth -$H$ during circulation, Eq. (26). The denominator contains an exponential term that accounts for thermal expansion along the wellbore (through $\mathrm{m}$, $\mathrm{n},$ and $\mathrm{p}$), while the numerator contains a root term that describes compression due to the increasing pressure (through $\mathrm{q}$, $\mathrm{r}$ and the local hydraulic resistance $\mathsf{\xi }(\mathrm{H})$). The coupling between temperature and pressure is intrinsically non-linear because $\mathsf{\xi }(\mathrm{H})$ itself depends on density and viscosity, both of which vary with $\mathrm{T}$ and $\mathrm{P}$.

For practical use, the implicit algebraic equations (29) and (32) must be solved. The laminar-flow model (29) is a polynomial of degree –q. Typically, $\mathrm{q}$ ≈ 6 for the aqueous phase of water-based fluids, see Table 5, because ${\mathrm{E}}_{\mathrm{up}}$ is of order ${10}^3$ MPa and its pressure derivative is of the order of 6. It can be solved iteratively for $\rho$ at a given depth and flow velocity. The accompanying inequality guarantees that the Reynolds number does not exceed the critical value 2300. Similarly, the turbulent-flow model (32) contains an additional term -$a\sqrt[4]{{\rho }^3}$ due to the Blasius friction factor; this term makes the equation slightly more non-linear but still manageable numerically.

3.3. Rheological Parameters Along the Wellbore

Once the density $\mathsf{\rho }(\mathrm{H})$ is known from either (29) or (32), the kinematic viscosity is obtained directly from the Frenkel-Andrade relation:

$$\mathsf{\nu }\left(\mathrm{H}\right)={\displaystyle \frac{\mathrm{A}\exp \left(\frac{\mathrm{E}}{\mathrm{k}\left({\mathrm{T}}_0+\mathrm{m}\mathrm{H}\right)}\right)}{\mathsf{\rho }\left(\mathrm{H}\right)}}.$$

(37)

The dynamic viscosity $\mathsf{\mu }(\mathrm{H})$ follows similarly. These profiles show that as temperature rises with depth, the exponential factor decreases strongly (viscosity drops), while the increase in pressure partly counteracts this drop by compressing the fluid and raising the viscosity. The net effect is a non-monotonic behaviour that can only be captured by the coupled model.

The shear stress in the annulus includes contributions from the hydrostatic pressure gradient, the frictional losses, and the local temperature and pressure gradients. For a typical deep well (e.g., $\mathrm{H}=4000$ m, $\mathrm{m}=0.025$ K/m, ${\mathrm{P}}_0=2$ MPa, $\mathrm{v}=1.2$ m/s), the model predicts that the shear stress at the wall increases by 15–20 % compared to an isothermal calculation, emphasizing the need for thermobaric correction.

3.4. Drilling Window and Accident Prevention

The integrated systems (35) and (36) provide a practical tool for drilling engineers. They simultaneously enforce:

The correct flow regime (laminar or turbulent one) via the velocity inequality;

The pressure balance in the annulus;

The condition that the bottomhole pressure stays between the formation pressure and the fracture pressure.

If the calculated ${\mathrm{P}}_{\mathrm{b}}$ from the density model falls outside this interval, the user can adjust the flow velocity $\mathrm{v}$ or the surface pressure ${\mathrm{P}}_0$ (or even change the fluid properties ${\mathsf{\rho }}_0$, $\mathrm{A}$, $\mathrm{E}$, etc.) until a safe operating point is found. The iterative solution is straightforward because the equations are monotonic in $\mathsf{\rho }$ and $\mathrm{v}$ over the range of interest.

3.5. Experimental Validation and Uncertainty

The experimental rig (Figure 1) was used to measure density and viscosity of water-based drilling fluids at temperatures from 1 °C to 110 °C and pressures from 10 MPa to 50 MPa. For a typical polymer-bentonite fluid, the measured densities at several depths were compared with the predictions of the laminar-flow model (29). The average relative deviation was 1.2 %, and the maximum deviation was 1.8 %, slightly above the claimed ±1.5 % but still acceptable for engineering purposes. The small discrepancy is attributed to uncertainties in the laboratory determination of $\mathsf{\alpha }({\mathrm{T}}_0,{\mathrm{P}}_0)$ and ${\mathrm{E}}_{\mathrm{up}}({\mathrm{P}}_0,{\mathrm{T}}_0)$, as well as to the assumption of a perfectly concentric annulus (eccentricity was not modelled).

The turbulent-flow model has not yet been validated against dedicated flow-loop or downhole measurements; however, it relies on the well-established Blasius friction law.

3.6. Order-of-Magnitude Assessment of the Competing Thermal and Baric Effects

The relative importance of the two competing effects can be assessed from Equation (26) without solving the full system. For a well of H = 4000 m with a temperature gradient of 0.0225 K/m, the temperature increases by 90 K along the annulus, so that the bottomhole temperature reaches 110 °C, which is exactly the upper limit of the range covered by the experiments. With a mean thermal expansion coefficient of approximately 5.2 × 10⁻⁴ K⁻¹, which corresponds to the measured range 20–110 °C (Table 1), the thermal term in the denominator of Equation (26) reduces the density by roughly 4.7%. The bottomhole pressure in such a well is of the order of 45–50 MPa, again within the range of the experiments; with a bulk modulus of about 2.2 × 10³ MPa (Table 3), the compression term increases the density by roughly 2%. The two effects are therefore of the same order of magnitude but of opposite sign, and the net change of density at the bottom is approximately −2.7%. For a fluid with a surface density of 1200 kg/m³, this corresponds to a decrease of about 32 kg/m³ at the bottom and, when integrated over the well depth, to a reduction of the hydrostatic pressure of the order of 0.5–1 MPa relative to the value calculated with the surface density. This estimate agrees in magnitude with the downhole density changes reported for water-based and synthetic-based fluids [5,17] and shows why neglecting either effect, or treating them independently, may lead to errors comparable with the width of a narrow drilling window.

3.7. Propagation of Input Uncertainties

Equation (26) also makes it possible to estimate how uncertainties in the input parameters propagate to the calculated density. Because ρ is directly proportional to ρ₀, the hydrometer error of up to 2 kg/m³ per 10 °C deviation from the reference temperature (Section 2.1) is transferred to the density profile without attenuation; for a fluid with a density of 1200 kg/m³, this corresponds to about 0.2%, which is small but not negligible compared with the thermobaric correction of several percent estimated in Section 3.6. The thermal term in the exponent of Equation (26) is proportional to the mean value of α over the temperature interval; hence, a relative error of 1% in the thermal expansion coefficient produces a relative density error of only about 0.05% under the conditions of Section 3.6. Likewise, a relative error of 1.5% in the bulk modulus leads to a density error of about 0.03%.

Consequently, the linearisation errors established in Section 3.1 (≤1% for α and ≤1.5% for E_up) contribute far less to the uncertainty of the predicted density than the accuracy of the surface density measurement and, above all, the uncertainty in the temperature profile along the annulus, which enters Equation (26) both linearly and quadratically through the gradient m. For example, a 10% error in the temperature gradient changes the calculated bottomhole density by about 0.5–0.6%, i.e., by 6–7 kg/m³ for the fluid considered above. This result emphasizes that reliable temperature data are the key prerequisite for the practical application of the model.

3.8. Worked Example for a Deep Well

To illustrate the practical use of the model, the calculation procedure of Section 2.9 was applied to a representative deep well. The input data are listed in Table 7. The depth interval and the temperature gradient were chosen so that the temperature and the pressure along the annulus remain within the ranges covered by the experiments described in Section 2.1 (up to 110 °C and up to 50 MPa): the bottomhole temperature of 110 °C and the bottomhole pressure of 48.7 MPa correspond to the upper limits of the tested ranges. The volumetric constants were obtained by least-squares fitting of the measured data: gradα(T) and α(T₀, P₀) from the row of Table 1 corresponding to 20 MPa, and grad(E_up) and E_up(P₀, T₀) from the row of Table 3 corresponding to 25 °C; the compressibility law is therefore extrapolated from 8 MPa to the pressures along the annulus, as discussed in Section 3.1 and Section 4.4. The viscosity constants correspond to a polymer–bentonite fluid with a dynamic viscosity of 20 mPa·s at surface conditions.

Because the Reynolds number exceeds 4000 over the whole depth interval, the turbulent branch, Equation (32), applies at every depth. The calculated profiles are shown in Figure 2, and the corresponding numerical values, which can also be used to reproduce the plots, are given in Table 8.

The density decreases monotonically from 1200 kg/m³ at the surface to 1165.4 kg/m³ at 4000 m, that is, by 2.9%, in close agreement with the order-of-magnitude estimate of Section 3.6 (−2.7%); the small difference arises because Equation (32) treats the coupling of the thermal, baric and frictional terms exactly, whereas the estimate evaluates them separately. The viscosity decreases by a factor of 3.6, so that the Reynolds number grows from 6.5 × 10³ to 2.3 × 10⁴ and the flow remains turbulent over the whole interval. The bottomhole pressure at 4000 m is 48.70 MPa, whereas a calculation with constant surface properties gives 50.46 MPa; the difference of about 1.8 MPa is equivalent to some 45 kg/m³ of mud weight and may be decisive when the drilling window is narrow.

The frictional contribution to the bottomhole pressure increases from 0.16 MPa at 500 m to 0.97 MPa at 4000 m and, in accordance with the sign convention adopted in Section 2.4, increases the bottomhole pressure, so that ECD* exceeds the local density of the fluid at every depth. It should also be noted that, with q ≈ 30, Equation (32) is a genuinely nonlinear algebraic equation; nevertheless, the iterative procedure of Section 2.9 converges to the adopted tolerance within a few iterations at every depth, which confirms that the model is suitable for real-time calculations.

4. Discussion

4.1. Comparison with Existing Models

Most earlier models [16,45,50] used separate multiplicative factors for temperature and pressure, i.e., $\mathsf{\mu }(\mathrm{P},\mathrm{T})={\mathsf{\mu }}_0\cdot \mathrm{f}(\mathrm{T})\cdot \mathrm{g}(\mathrm{P})$. Such an approach assumes that the two effects are independent, which is not strictly true because the compressibility $\mathsf{\beta }$ itself depends on temperature and the thermal expansion $\mathsf{\alpha }$ depends on pressure. The present model inherently couples $\mathrm{T}$ and $\mathrm{P}$ through the integration of the total differential (6) and through the pressure-dependent –$\mathsf{\xi }(\mathrm{H})$ that contains viscosity. In addition, the Herschel–Bulkley parameters are indirectly incorporated via the choice of the laminar or turbulent friction factor. This simplification is adequate for many engineering applications, but it should be examined in future works.

A notable advantage of the proposed model is that it yields closed-form (albeit implicit in ρ) algebraic equations, avoiding the need for time-consuming numerical integration of differential equations along the wellbore. This makes it suitable for real-time drilling optimisation.

4.2. Physics-Based and Data-Driven Density Models

Data-driven models of downhole density achieve very low errors on their test datasets [24,25,26,27], but they require large and representative training sets and act as black boxes with respect to the underlying physics. The proposed model occupies the opposite end of the spectrum: it is derived from the thermodynamic definitions of the thermal expansion and compressibility coefficients, contains only five fluid-specific volumetric constants (ρ₀, gradα, α(T₀, P₀), grad(E_up) and E_up(P₀, T₀)) and two viscosity constants, and each constant can be measured with standard laboratory equipment. Its accuracy is limited mainly by the validity of the linear approximations, which was verified statistically in Section 3.1.

The two approaches are complementary rather than competing. On the one hand, the closed-form model can generate physically consistent synthetic data and provide bounds for the output of machine-learning models, reducing the risk of non-physical predictions outside the training range. On the other hand, data-driven correction terms trained on field measurements may compensate for effects not included in the model, such as dissolved gas, variable solids content or thermal transients after circulation breaks. Machine-learning methods have already been applied successfully to related drilling problems, including real-time estimation of rheology and wellbore stability assessment [14,62], and their combination with thermobaric density modelling appears promising.

Compared with compositional models [17,19], which require separate volumetric data for each liquid component and the volume fractions of all phases, the proposed approach characterises the fluid as a whole by two linear laws. This makes it less general but considerably simpler to parameterise for a specific fluid system, which is advantageous in field applications where the composition changes during drilling and the parameters must be updated quickly. Automated laboratory measurements of the density of drilling fluids as a function of pressure and temperature [20] provide exactly the type of data required for such rapid re-parameterisation.

4.3. Integration with Real-Time Drilling-Fluid Monitoring

The computational simplicity of the model makes it suitable for integration into automated drilling-fluid monitoring systems. Modern surface measurement units provide density and rheological readings at intervals of minutes [62,63,66], and software–hardware complexes developed for continuous monitoring of drilling-fluid rheology can transmit these data to the drilling-control computer [32,67]. In such an architecture, the model acts as a physics-based layer that transforms each new surface measurement of density and viscosity into the downhole density and pressure profiles along the annulus. Because the equations are solved in milliseconds, the transformation can be repeated after every measurement cycle, and the resulting ECD can be compared with the drilling window in real time.

This capability is particularly valuable in MPD operations, in which the surface back-pressure P₀ is actively controlled to keep the bottomhole pressure within narrow limits [54,64,65]. Because P₀ appears explicitly in Systems (35) and (36), the model can be used not only to check the current state but also to calculate the back-pressure required to compensate for thermobaric density changes when the circulation rate is changed. The case of stopped pumps, however, requires an extension of the model to non-circulating conditions with a transient temperature field, as discussed below.

More broadly, the model fits into the concept of digital twins of the drilling-fluid system, in which surface measurements, hydraulic models and downhole data are continuously reconciled [69]. Within such frameworks, the closed-form thermobaric model can serve as a fast forward model for the density component, whereas systematic deviations of its predictions from downhole pressure-while-drilling data can indicate anomalies such as influxes, losses or changes in fluid composition [35].

4.4. Limitations and Future Work

Several limitations of the present model should be acknowledged. First, the model assumes a concentric annulus, whereas eccentricity can significantly alter frictional pressure losses and is not accounted for. Second, the Frenkel-Andrade equation for viscosity is valid only for Newtonian or slightly non-Newtonian fluids; strongly shear-thinning fluids would require a more detailed rheological model, such as a Herschel-Bulkley model coupled with the thermobaric dependencies. Third, the linear approximations for α(T) and E_up(P) were verified only for water-based fluids; for oil-based or synthetic fluids, the linearity ranges may differ. Finally, the model does not consider chemical reactions or phase changes (e.g., gas evolution) that could occur under extreme HP/HT conditions. Future research should therefore extend the model to handle eccentric annuli and non-Newtonian yield-stress fluids, and should validate the turbulent flow equations with dedicated experiments on a full-scale flow loop.

In addition, the model uses a constant temperature gradient, whereas the actual temperature of the circulating fluid differs from the undisturbed geothermal profile and depends on the circulation time, the flow rate and the heat exchange between the drill string and the annulus; measurements with drilling microchips and other downhole sensors show that these differences can be substantial [4]. The transitional flow regime (2300 < Re < 4000) is treated only approximately, and the effects of drill-string rotation and eccentricity on frictional losses, which are known to be significant [6,61], are neglected. The linear law E_up(P) was established from measurements at pressures up to 8 MPa and is extrapolated to the 10–50 MPa range in all calculations; the measured slope, grad(E_up) ≈ 12–30, is also larger than the pressure derivative of the bulk modulus of pure water reported in the literature, which is of the order of 5–7 [22]. Repeating the volumetric tests over the full pressure range of the rig is therefore the first priority for further work, since the value of this slope directly controls the compression term of Equation (26). Finally, the linearity of α(T) and E_up(P) was established for water and should be verified for fluids with high solids content, high salinity or significant amounts of dissolved gas. These points define the main directions for further development of the model.

5. Conclusions

In this study, a comprehensive thermobaric model of drilling fluid in the annular space of deep wells has been developed. Using the total differential of a function of two variables together with experimental data for wellbore water, it was first established that, at pressures of 10 MPa and higher, the coefficient of thermal expansion $\alpha$ is a linear function of temperature with an error below 1 %, and the isothermal bulk modulus ${\mathrm{E}}_{\mathrm{up}}$ is a linear function of pressure with an error not exceeding 1.5 % at temperatures of 5 °C and above, within the investigated pressure range of 0.5–8 MPa; the latter law is extrapolated to the pressures of deep wells. These linear approximations, validated through correlation coefficients, allowed the derivation of a differential equation describing the change in drilling fluid volume under simultaneous thermal and baric loading.

The differential equation was integrated to obtain a closed-form expression for the fluid density $\mathsf{\rho }$ at any depth $\mathrm{H}$ during circulation. The resulting density formula explicitly accounts for the non-linear interaction between temperature and pressure through the combined effects of thermal expansion, pressure compression, and hydraulic resistance in the annulus. Two practical variants of the density equation were derived: one for laminar flow and one for turbulent flow. Each variant reduces to an algebraic equation coupled with a Reynolds-number inequality that guarantees the prescribed flow regime. The model also provides depth profiles of kinematic viscosity and shear stress, using the Frenkel-Andrade relation for dynamic viscosity.

To ensure safe drilling operations, the density model was embedded within the drilling window constraints. The final systems of equations and inequalities (one for laminar, one for turbulent flow) allow the engineer to select the circulation velocity and fluid parameters that keep the bottomhole pressure between the formation pressure and the fracture pressure, thereby preventing lost circulation, well collapse, and pipe sticking. The models are intended for wells deeper than 1000 m, where the temperature exceeds 20 °C and the pressure exceeds 10 MPa, and their overall error does not exceed ±1.5 % provided that the input parameters are determined with the accuracy stated in Section 3.7 and that the linear law E_up(P), established here up to 8 MPa, remains valid at higher pressures.

The practical value of the proposed thermobaric model lies in its ability to provide rapid and physically consistent predictions of downhole density and pressure for the hydraulic design of deep and ultra-deep wells under HP/HT conditions. Owing to its low computational cost and small number of measurable parameters, it can be integrated into real-time drilling-fluid monitoring and MPD control systems. Further work will focus on extending the model to non-circulating conditions with a transient temperature field, to eccentric annuli and yield-stress fluids, and to oil-based and synthetic drilling fluids.

6. Patents

1. Muratova, S.K.; Ratov, B.T.; Khomenko, V.L.; Kamyshatskyi, A.F.; Seidaliev, A.A.; Zhailiev, A.O.; Omirzakova, E.Zh.; Aliakbar, M.M. System for Automatic Measurement of Drilling Mud Density. KZ Patent 37868, 6 March 2026.

2. Ratov, B.T.; Khomenko, V.L.; Kamyshatsky, A.F.; Muratova, S.K.; Syrlybekkyzy, S.; Zhailiev, A.O.; Omirzakova, E.Zh.; Makyzhanova, A.T. Device for Automatic Measurement of Rheological Properties of Drilling Mud. KZ Patent 37993, 17 April 2026.