Automatic Measurement of the Dependence on Pressure and Temperature of the Mass Density of Drilling Fluids

Drilling fluid mass density is essential for the calculation of the hydrostatic pressure in a borehole. It is also an important parameter in viscous friction calculations because it appears in the Navier-Stokes equation and in buoyancy calculations when performing torque and drag or drillstring dynamic analysis. During drilling operations, the drilling fluid mass density is usually measured with a mud balance at regular intervals. The precision of the mud balance measurements is about ±10 kg/m³. Furthermore, drilling fluids are subject to a wide range of pressures and temperatures, and therefore the drilling fluid mass density cannot be assumed to be constant for the range of depths of normal oilfield wells. Isambourg et al. (1996) found that the pressure and temperature dependence of drilling fluids can be described by a biquadratic function in the temperature range of 20–200°C and a range of pressure between 1 and 1400 bar. The biquadratic model proposed by Isambourg has since been used in a more recent study (Hermoso et al. 2017). Based on multiple published density measurements made on various fluids (Zamora et al. 2000, 2012; McMordie et al. 1982; Peters et al. 1990; Hemphill and Isambourg 2005; Demirdal and Cunha 2007, 2009; Demirdal et al. 2007; Hussein and Amin 2010), it has been found that the PVT model can be reduced to a linear function of temperature combined with a quadratic form of pressure without losing much of its accuracy and which is the model published in API RP 13D (2017):

where $ρm$ is the mass density, $T$ is the temperature, $p$ is the pressure, and $Am$⁠, $Bm$⁠, $Cm$⁠, $Dm$⁠, $Em$⁠, and $Fm$ are the parameters of the PVT model. The accuracy of the model described by Eq. 1 is about 0.35% on average with a standard deviation of 0.23% for a temperature range from 3 to 260°C and pressures from 1 to 2000 bar, based on measurements of 26 different fluids.

Calibrating Eq. 1 for a particular drilling fluid requires access to a high-pressure and -temperature densitometer, which is not routinely available at the rigsite, except, perhaps, during drilling operations in high-pressure, high-temperature conditions. The alternative is to estimate the parameters of the PVT model by combining the PVT behavior of each component of the drilling fluid. Because a drilling fluid is a mixture of solids and liquids, its mass density is simply the volumetric weighted average of all the components:

where $ρm$ is the drilling fluid mass density, $Ω$ is the set of indices for the different components, $fi$ is the volume fraction of the ith component, and $ρi$ is its corresponding mass density.

The PVT behavior of several aqueous solutions of NaCl, NaBr, KCl, KBr, CaCl₂, CaBr₂, ZnCl₂, and ZnBr₂ has been described by Kemp et al. (1989). Their model is relatively complex and does not reduce to Eq. 1. For that reason, Eq. 1 has been expanded to account for the concentration of salt by replacing the parameter $Ao$ with a cubic form of the weight of salt:

where $ρw$ is the mass density of the brine and $wt$ is the weight fraction of the salt at the normal condition of temperature and pressure (atmospheric pressure and 297 K). $S0$⁠, $S1$⁠, $S2$⁠, $S3$⁠, $Bw$⁠, $Cw$⁠, $Dw$⁠, $Ew$⁠, and $Fw$ are the coefficients of the model as published in API RP 13D (2017).

The PVT behavior of base oils can also be described by Eq. 1. Zamora et al. (2013) presented parameters for several base oils and some brines. Moreover, the solid components have coefficients of thermal expansion and compressibility (Isambourg et al. 1996) which should be accounted for. Knowing the PVT characteristics of each component, it is possible to estimate the mass density of the fluid mixture with Eq. 2. However, as the various components have different compressibilities and thermal expansion characteristics, the volumetric fraction of each component changes with pressure and temperature. A method to calculate changes in volume fractions with pressure and temperature is described by Cayeux et al. (2014).

In practice, the downhole drilling fluid density will also be influenced by the presence of cuttings, dissolution of formation solids, absorption of formation fluids in the liquid phase, and barite sagging. These additional considerations are not addressed in this paper because the objective is to propose an alternative system of surface measurement to those which are, themselves, subject to the same limitations.

However, the problem remains that the composition of the drilling fluid changes continuously throughout a drilling operation, either because of the actions taken by the fluid engineer to maintain the quality of the drilling fluid or because formation fluids and solids are entrained during circulation. Even more problematically, it is often not easy to get an accurate description of the component quantities of the drilling fluid. Cayeux (2020) illustrated the impact of this lack of information on the estimated hydrostatic pressure. For instance, if the drilling fluid density is known at atmospheric pressure and 50°C (e.g., 1700 kg/m³) and nothing else, the 95 percentile estimate of the hydrostatic pressure at 3500 mTVD (true vertical depth), lies between −5 and +7 bar (Fig. 1a). However, for an oil-based fluid (OBM) with a 70/30 oil/water ratio (OWR), the 95 percentile at 3500 mTVD reduces to −2.5 and +3 bar (Fig. 1b). Knowing the fluid density for a water-based drilling fluid (WBM) at three different temperatures under atmospheric conditions will reduce the confidence interval, but only to −4 to 5.5 bar (Fig. 1c).

Having to rely on manual measurements and assumptions about the drilling fluid formulation presents difficulties when trying to automate the drilling process. Therefore, it would be a considerable advantage if the drilling fluid mass density could be measured in real time. Although such information is available on some rigs, it would be even better to obtain a full description of the drilling fluid PVT characteristics directly from inline measurements. Therefore, the goal of the present work is to establish whether it is possible to design an inline apparatus capable of measuring the PVT behavior of a drilling fluid with acceptable precision (for example, to ±5%) over a wide range of pressures and temperatures, typically 10–200°C and from 1 to 2000 bar.

Because Eq. 1 has six parameters, it is necessary to have at least six measurements made under different conditions of pressure and temperature. An obvious approach would be to measure the mass density of the fluid under six different pairs of pressure and temperature; consequently, from Eq. 1, the following set of simultaneous equations can be written:

where $ρm1$⁠, $ρm2$⁠, $ρm3$⁠, $ρm4$⁠, $ρm5$⁠, and $ρm6$ are the measured mass densities at the respective conditions of pressure and temperatures $p1,T1$⁠, $p2,T2$⁠, $p3,T3$⁠, $p4,T4$⁠, $p5,T5$⁠, and $p6,T6$⁠. This system of equations is linear in $Am$⁠, $Bm$⁠, $Cm$⁠, $Dm$⁠, $Em$⁠, and $Fm$ and can be easily solved using a Gauss-Jordan method, for instance (Teukolsky et al. 1992).

It is desirable to avoid operating at high pressures and temperatures to keep the design of the apparatus as simple as possible. Considering that the model of Eq. 1 has an average accuracy of about 0.35% over a range of temperatures spanning more than 250°C and a range of pressures up to 2000 bar, it is possible to calibrate the parameters with measurements made on a smaller range of pressures and temperatures. However, the accuracy of the calibrated model will depend on the accuracy of measurements and the range of pressures and temperatures over which they are taken. Fig. 2 illustrates how measurement uncertainty and the associated range of measurements influence the uncertainty of the extrapolated value when utilizing a model of known precision. In this illustration, the model is calibrated using randomly chosen values of the measurements within their accuracy range to provide a stochastic description of the envelope of model predictions outside the measurement range. The precision of the measurements limits the acceptable extrapolation range, and a narrow measurement range results in increasing uncertainties the greater the extrapolation.

Accounting for model accuracy and measurement precision and repeatability, it is useful to establish a calibration range of pressures and temperatures for the sample measurements that will permit extrapolation of the model to higher pressures and temperatures with acceptable accuracy. If the densitometer has a measurement precision of 0.05 kg/m³ and repeatability of 0.01 kg/m³, based, for example, on an Anton Paar L-dens 7500 inline densitometer, stochastic simulations can be run to quantify the measurement error for a suitable calibration range of pressures and temperatures. Measurement precision is regarded as a source of systematic error (i.e., one randomly drawn value is applied to all measurements), while measurement repeatability is considered as a random error drawn individually for each measurement. For example, if the PVT model is to be used over a range of 1000 bar and 200°C, the root mean square (RMS) of the difference between the predicted density and the “true” value can be calculated with random variations in the systematic bias of the density measurement when the calibration samples lie within a small range of pressure and temperature. The results of 250 stochastic simulations presented in Fig. 3 show that to achieve a 5% average accuracy for the entire range of 1000 bar and 200°C, it is necessary to use a calibration range of at least 30 bar and 25°C.

A possible design for an inline apparatus could involve pumping the drilling fluid past a controllable heating element with a controllable choke downstream of a densitometer applying backpressure while measuring the mass density (Fig. 4). The set points for the heating element and the choke could be changed six times to collect the necessary mass densities to calibrate the PVT model of Eq. 1. However, in practice, the measurement sequence would take too long because changing the heating element temperature would require time to equilibrate. Gathering a complete set of calibration measurements might easily take 15–30 minutes. In practice, with a typical annular velocity of 0.8 m/s, if PVT measurements are available every 30 minutes, then the same measured values would apply for an annular fluid column of about 1400 m in length. In a typical deep well, it might mean only four or five PVT measurements in the entire annulus. This would affect the resolution of the bottomhole pressure estimate, especially if the device had measured a momentarily heavy or light sample from the mud system. Conversely, if the PVT measurements are available every 30–60 seconds, with typical annular velocities around 0.8–1 m/s, the same PVT measurement will apply for about 30 m (i.e., one stand), which is a reasonable discretization length for hydraulic calculations. So, the aim of the inline PVT measurement device should be to provide measurements at most every minute.

To speed up measurements, an alternative would be to perform six measurements simultaneously. Fig. 5 shows a configuration where the densitometers are mounted in series. The configuration would consist of six parallel branches or any combination of series and parallel branches. Fig. 6 depicts two parallel branches where, in one branch, the temperature of the fluid remains constant, while in the second, it is changed. For each of the branches, backpressure is applied at two intermediate positions. This configuration has the advantage of fewer chokes and pressure sensors, four instead of five.

With simultaneous measurements of mass density, it is possible to provide in real time a complete, calibrated, PVT model. Yet, the circuitry is complex, involving four or five chokes. Therefore, it would be worthwhile to investigate whether the hydraulic circuit could be further simplified.

The pressure and temperature dependence of drilling fluid densities is related to the fluid compressibility and its thermal dilatation. Therefore, if it is possible to directly measure the fluid compressibility, it would be possible to simplify the setup and reduce the number of pressure and temperature conditions required for the solution. The isothermal compressibility (⁠$κT$⁠) of a fluid is defined as follows (Landau and Lifshitz 1980):

where $V$ is the volume that is related to the mass density through

where $m$ is the mass of matter contained in the volume $V$⁠. It is possible to express the isothermal compressibility of the fluid in terms of the PVT model described in Eq. 1:

Let us now consider the definition of the speed of sound (⁠$c$⁠) in matter:

where $S$ is the entropy. Then, the isentropic compressibility is defined as follows:

and it is possible to express a relation between the isentropic compressibility to the speed of sound in the medium:

Finally, the isentropic compressibility is related to the isothermal compressibility through the following relation (Landau and Lifshitz 1980):

where $α=1V(∂V∂T)p$ is the volumetric coefficient of thermal expansion and $Cp$ is the isobaric specific heat capacity. It is possible to express the volumetric coefficient of thermal expansion in terms of the PVT model described in Eq. 1:

Thus, if the mass density and the speed of sound are measured at the same time, six equations can be obtained with only three combinations of pressure and temperature:

where $c1$⁠, $c2$⁠, and $c3$ are the respective speeds of sound for the conditions of pressure and temperature $p1,T1$⁠, $p2,T2$⁠, and $p3,T3$⁠, $κp1,T1$⁠, $κp2,T2$⁠, and $κp3,T3$ are the estimated isothermal compressibility coefficients at the respective conditions of pressure and temperature $p1,T1$⁠, $p2,T2$⁠, and $p3,T3$⁠, and $αp1,T1$⁠, $αp2,T2$⁠, and $αp3,T3$ are the estimated volumetric coefficients of thermal expansion at the respective conditions of pressure and temperature $p1,T1$⁠, $p2,T2$⁠, $p3,T3$⁠.

The resulting system of equations is nonlinear. To solve that system of equations, Eq. 13 is reformulated:

where $x=(AmBmCmDmEmFm)$ , $Xj=(0000000Xm0000000000000Xm0000000000000Xm)$ with $j=1$⁠, $Xm=Bm$⁠, $j=2$⁠, $Xm=Dm$⁠, $j=3$⁠, $Xm=Fm$⁠. The $6×6$ matrix $A$ is as follows:

the matrix $B$ is

the matrix $C$ is

the matrix $D$ is

and $Et=(1c121c221c32ρm1ρm2ρm3).$

Eq. 14 can then be rearranged as follows:

Noting that $X1$⁠, $X2$⁠, and $X3$ are expressed as a function of $x$⁠, Eq. 19 has the form $x=fx$⁠, which can be solved efficiently by a fixed-point iteration method:

where $j$ denotes the iteration level. The algorithm stops when $xj+1-xj≤ε$⁠, $ε$ being a small real value.

Because the term $α2TρCp$ in Eq. 11 is much smaller than $κT$⁠, it is therefore possible to use

as an initial value for the fixed-point iteration method. In practice, the solution of Eq. 13 with the fixed-point iteration method does not require more than a half dozen iterations to achieve a stopping condition of $ε=10-6$⁠.

Suppose the speed of sound can be measured with an uncertainty of 0.1 m/s and repeatability of 0.01 m/s, based, for example, on the characteristics of an inline combined densitometer and speed-of-sound sensor model L-Com 5500 from Anton Paar, then it is possible to proceed as before and estimate the range of pressures and temperatures required to obtain a predictability of the mass density over a range of 1000 bar and 200°C. The results in Fig. 7 show that it is possible to achieve a precision better than 5% on average over the whole range of pressure and temperature of 1000 bar and 200°C, with three calibration measurements over a range of 60 bar and 5°C.

This result shows that it is possible to reduce the number of pressure and temperature conditions to three by measuring the speed of sound in the drilling fluid in addition to the mass density. However, this reduction in the number of pressure and temperature conditions depends on knowing the specific heat of the drilling fluid, which is often not the case. In the next section, the influence of the fluid specific heat capacity on the heat transfer along an insulated pipe will be analyzed.

The definition of specific heat is (Landau and Lifshitz 1980) as follows:

where $Q$ is the amount of heat (i.e., it characterizes the amount of energy necessary to change the temperature of a given mass of matter).

The specific heat capacity of the drilling fluid can be obtained using a hydraulic circuit consisting of a pump, a volumetric flow-rate measurement, and an insulated pipe in which a drilling fluid is circulated with a heating element upstream of the pipe (Fig. 8).

Consider the energy conservation of a system under isobaric conditions (Corre et al. 1984):

where $H$ is the enthalpy per unit mass, $Q→f=ρmHv→$ is the forced-convective heat flow, $Q→c$ is the conductive and natural convective heat flow, and $qs$ is a heat source per unit time and per unit length.

For an isotropic solid, like steel or an insulating material, the pure conductive heat flow is expressed as follows:

where $λ$ is the thermal conductivity.

The Nusselt number in the heat flow direction (⁠$NuL$⁠) is introduced to estimate the natural convective heat flow in a drilling fluid:

where $hL=qΔT$ is the convective heat transfer coefficient in the heat flow direction, $L$ is the characteristic length, $q$ is the heat flux (i.e., thermal power per unit area), and $ΔT$ is the temperature difference. For laminar flow, Corre et al. (1984) use the following expression of the Nusselt number (Corre et al. 1984):

where $n1=πλmL4QvρmCpm$⁠, $Qv$ is the volumetric flow rate, $Cpm$ is the specific heat capacity of the drilling fluid, $λm$ is the thermal conductivity of the drilling fluid, and $ψ$ is defined as follows:

Furthermore, the convective wall heat flux in laminar flow conditions (⁠$qw$⁠) is then (Barletta 1996)

where $Ri$ is the internal pipe radius, $Tw$ is the wall temperature, and $Tm$ is the drilling fluid temperature.

The enthalpy of a thermodynamic system is the sum of its internal energy $U$ and the work required to achieve its pressure and volume:

In differential form, Eq. 29 can then be expressed as [see Landau and Lifshitz (1980) for a demonstration]:

Assuming that the pressure drop along the pipe length of Fig. 8 is negligible, then $dH≈CpdT$⁠, so the energy balance of the drilling fluid contained in the pipe can be expressed as:

where $Tm$ is the drilling fluid temperature, s is the curvilinear abscissa, and $Ti$ is the temperature of the wall inside the pipe. The characteristic length $L$ of the Nusselt number is the internal pipe diameter (i.e., $L=2Ri$⁠). The hydraulic work arising from viscous pressure losses is

The energy balance for the pipe and the insulation material is defined in a cylindrical coordinate system as

where $ρs$⁠, $Cps$⁠, and $λs$ are, respectively, the mass density, specific heat capacity, and thermal conductivity of the solid material of either the pipe or the insulation material, and $qsh$ is the generated heat per unit time and unit length of the heating element. Considering the axial symmetry of the system, the term $∂2Ts∂r2$ is necessarily equal to zero.

In steady-state conditions, without heat generation and in pure axial conduction, Eq. 33 can be written as follows:

After integration over a cylindrical area between $R1$ and $R2$⁠, the axial heat flux per unit length (⁠$qa$⁠) is derived:

The axial heat flux over an axial distance $l$ can be estimated:

where $∆Ta$ is the axial difference of temperature.

Similarly, in steady-state conditions without heat generation and in pure radial conduction, Eq. 33 can be written as follows:

After integration over a shell of radius $r$ and length $l$⁠, the radial heat flux per unit radial distance (⁠$qr$⁠) is as follows:

And therefore, the radial heat flux between $R1$ and $R2$ is

where $∆Tr$ is the radial temperature difference.

Analysis of the heat transfer along an insulated pipe after the fluid has passed a heat source shows that the temperature along the pipe depends not only on the specific heat capacity but also the thermal conductivity of the fluid.

Attaching temperature sensors to the outside of the pipe in Fig. 8 enables the radial heat flux to be measured while circulating fluid in the system. This proposition was verified by constructing a flow loop similar to that shown in Fig. 8. Fig. 9 shows a schematic view of the experimental setup. A 10-L tank is equipped with an agitator and copper coil through which tap water is circulated. The flow rate is controlled so that the temperature in the tank remains constant.

A progressive cavity pump injects the fluid into the insulated instrumented pipe section. The pump is equipped with an encoder to measure the pump speed, and because there is very little pressure loss in the flow loop, the flow rate is directly proportional to the pump speed. The conversion factor from pump speed to volumetric flow rate has been calibrated by measuring the time to pump a known volume of fluid. The insulated pipe section comprises a 2-m length of 20-mm-diameter stainless-steel tube. Close to the inlet is a cylindrical aluminum block 200 mm in length and 100 mm in diameter heated by four resistors. A temperature sensor (TT1) is placed inside the aluminum block, and a controller regulates the electrical current to the heating resistor to regulate the temperature of the aluminum to a given set point. The remaining temperature sensors are placed 250 mm apart, and the insulation material is 30 mm thick.

Fig. 10 shows what happens when the heating element is turned on at 550 seconds, and then just after 3,100 seconds, when the flow rate is decreased from 0.451 to 0.096 L/min. The first temperature sensor downstream of the heating element (TT2 in Fig. 10) is positioned 100 mm above the top of the heating element and therefore is dominated by the heat conducted through the pipe with minimal influence from the circulating fluid. Because the mass density is known, temperature measurements at two different fluid velocities allow the specific heat capacity and the thermal conductivity of the drilling fluid (Fig. 11) to be derived from the heat transfer model described by Eqs. 31 and 33.

The measurement of radial heat flux along a pipe under varying flow conditions, shown in Fig. 12, with three stages of pressure and temperature permits the determination of the PVT properties of the fluid, as follows:

1. Ambient pressure and temperature conditions at the inlet

2. Increasing the temperature by 5–10°C at a pressure in the region of 50–60 bar

3. Increasing the temperature by a further 5–10°C with a pressure in the region of 25–30 bar

The volumetric flow rate and mass density at the inlet are measured with a Coriolis flowmeter. A speed-of-sound sensor is installed upstream of the inlet of the progressive cavity pump. In the second stage, a heating element downstream of the pump increases the fluid temperature by 5–10°C. A stainless-steel pipe is mounted on a thermal insulator downstream from the heating element. A differential pressure sensor measures the pressure losses along the pipe. The temperature distribution along a pipe is measured by at least four temperature sensors after which the mass density and speed of sound are again measured with a Coriolis flowmeter. The pressure in this second stage is controlled by a choke downstream of the flowmeter. The third stage consists of a pipe of a different diameter from the first, designed to create flow conditions that differ from those in the previous stage. A similar heating element further increases the temperature. The pipe is thermally insulated from the heating element and, again, at least four temperature sensors are placed along the pipe to measure the temperature distribution as well as a differential pressure transducer. Another combined sensor measures the mass density and the speed of sound under the third set of conditions. The pressure in this third stage is controlled by a choke at the outlet.

The heat transfer (Eq. 31) depends on $ρm,$$Cpm$⁠, $λm$⁠, and the Nusselt number $NuL$⁠, which in turn depends on $ρm,$$Cpm$⁠, $λm$⁠. The solution of Eq. 31 gives the temperature evolution of the system. Because the fluid’s mass density is measured, it remains to determine $Cpm$ and $λm$⁠. These parameters are obtained by solving the minimization problem:

where $Ti$ and $T~l$ are, respectively, the measured and estimated temperatures at $nT$ different positions along the pipe. If Eq. 31 is divided by $Cpm$⁠, then the equation depends on a term proportional to $λmNuLCpm$ and $qsΔpCpm$⁠. Because the second term is small compared with the first, the temperature along the pipe is mostly influenced by the ratio of $λm$ and $Cpm$⁠. Because it is desired to explicitly determine the specific heat capacity of the fluid, it is necessary to create conditions that yield both the specific heat capacity and the thermal conductivity. Two dissimilar pipe sizes will result in different fluid velocities, allowing the term $qsΔpCpm$ to elicit this difference.

To extract values of C_p and λ from the setup shown in Fig. 12, the system is discretized as shown in Fig. 13 in $na$ axial elements of length $∆s$ (indexed by $i$⁠) and in radial elements of radial distance $∆r$ (indexed by $j$⁠). Similarly, time may be discretized in time steps of duration $∆t$ (indexed by $k$⁠).

Eq. 31 governs the heat transfer in the drilling fluid and in the air (outside of the system). For air, there are no heat sources, and there is no forced circulation; therefore, the energy balance at the external wall of the setup reduces to

where $ρair$ is the mass density of air, $Cpair$ is the specific heat capacity of air, $Rnr$ is the external radius, $hair$ is the air convective heat transfer coefficient, and $Tair$ is the air temperature.

For the drilling fluid, heat generation arises from viscous pressure drops. Eq. 32 is used to estimate the heat generation caused by the passage of the fluid through the chokes. The energy balance at the internal wall of the pipe is

where $R1$ is the internal radius of the pipe, $Tmi$ is the temperature of the fluid at axial position i, and $∆pi$ is the viscous pressure drop at axial position i. Note that $Qv$ and $NuL$ also depend on the axial position.

For any other values of $j$⁠, the energy conservation is based on Eq. 33:

where $Rj+12$ means the average radius between $j$ and $j+1$ and $qshi$ is the heat generation at axial position i, if any.

After regrouping all the terms indexed $k+1$ on one side and the terms indexd $k$ on the other, a system of linear equations is obtained. It can be solved for all $Ti,j,k+1$⁠, for example, using a Gauss-Jordan method:

where $P$ is a pentagonal matrix in which coefficients depend on the current boundary conditions and values estimated at the previous timestep $k$⁠, $Tk+1$ is a vector of unknown temperature at timestep $k+1$⁠, and $R$ is a vector of coefficients that depend only on the current boundary conditions and values estimated at the previous timestep $k$⁠. The boundary conditions are the air temperature, the temperature of the pipe and the drilling fluid at the inlet, the volumetric flow rate at the inlet, the drilling fluid mass density in in-situ conditions, the specific heat capacity and the thermal conductivity of the drilling fluid, and the heat generated by the heating elements. In practice, a Crank-Nicolson scheme is used to improve the stability of the numerical resolution (Crank and Nicolson 1947). It is then possible to calculate, as a function of time, the temperature anywhere in the system and on the external side of the pipe.

Initial values of $Cpm$ and $λm$ are unknown and therefore must be estimated. Because $nT$ temperatures are measured along the pipes, $Cpm$ and $λm$ can be estimated by solving Eq. 40.

Setting $F=∑l=1nTTl-T~l2=∑l=1nTrl2$⁠, the sum of squared temperature differences, and defining $βt=[Cpmλm]$⁠, the vector of model parameters. $F$ is a minimum when

Linearizing the solution vector from one iteration to another:

where $ξ$ is the iteration number. Then the first order approximation of $T~lβ$ is

where $Jls=∂rl∂βs$ are the members of the Jacobian matrix $J$⁠. Then, a first-order approximation of $rl$ is

where $∆Tl,ξ=Tl-T~lβξ$ is the temperature difference for sensor $l$ at iteration $ξ$⁠. Substituting in Eq. 45, here is the resulting equation: