Hybrid Newton–Successive Substitution Method for Multiphase Rachford-Rice Equations

In multiphase (≥3) equilibrium calculations, when the Newton method is used to solve the material balance (Rachford-Rice) equations, poorly conditioned Jacobian can lead to false convergence. We present a robust successive substitution method that solves the multiphase Rachford-Rice equations sequentially using the method of bi-section while considering the monotonicity of the equations and the locations of singular hyperplanes. Although this method is slower than Newton solution, as it does not rely on Jacobians that can become poorly conditioned, it can be inserted into Newton iterations upon the detection of a poorly conditioned Jacobian. Testing shows that embedded successive substitution steps effectively improved the robustness. The benefit of the Newton method in the speed of convergence is maintained.


Introduction
Petroleum reservoir fluids are multi-component mixtures primarily made of hydrocarbons, held subsurface at high-temperature and high-pressure conditions [1,2]. Although under most situations the phase behavior of these fluids can be adequately described by vapor-liquid two-phase equilibria, multiphase phenomena (number of phases ≥3) have also been widely observed. Some of them, such as wax precipitation [3][4][5][6], hydrate formation [7][8][9][10], and asphaltene precipitation [11][12][13][14] are of significant importance to oil and gas production and transportation. In enhanced oil recovery, mixtures of reservoir oil and injected gas, such as CO 2 , can also exhibit complex phase behaviors at low pressures [15][16][17]. For these reasons, research on multiphase equilibrium calculations is still very active.
Multiphase (≥3) equilibrium calculation includes phase-stability [18] and phase-split [19] calculation steps. These two steps are usually carried out in series, except for the class of methods that was started by Gupta et al. [20], where they are carried out simultaneously. In the phase-split calculation step, mass balance equations and fugacity equations are solved. These two sets of equations can be solved sequentially using the method of successive substitution: the mass balance equations are first solved for a given set of equilibrium ratios; then, the equilibrium compositions of the phases that were obtained from the solution of the mass balance equations are used to update the equilibrium ratios through each phase's equation of state or activity correlation [21][22][23]. These two sets of equations can also be solved simultaneously by the Newton method [24]. Although the Newton method converges more rapidly than the method of successive substitution, its accuracy is often influenced by the initial guess. For this reason, Michelson [19] recommended to begin phase-split calculations with successive substitutions and then switch to the Newton method. Successive substitution method and Newton method for phase-split calculations are abbreviated as SS-PS and NM-PS in this study.

Successive Substitution Method for Multiphase Rachford-Rice Equations
In this study, we use n j to denote the mole fraction of phase j in a multiphase, multicomponent mixture, x j i to denote the mole fraction of component i in phase j, and z i to denote the mole fraction of component i in the mixture. N C is the total number of components and N P is the total number of phases. The equilibrium ratio of component i between phase α and phase β is defined as K αβ i = x α i /x β i . The multiphase Rachford-Rice equations that are derived from mass balances are Here, selection of reference phase "1" is arbitrary, but it should contain every component of the system. While the selection of reference phase does not affect the results, it does affect convergence path and speed. We recommend to select reference phase to avoid extreme values of equilibrium ratios. For a given set of equilibrium ratios K j1 i (i = 1, 2, · · · , N C , j = 2, · · · , N P ) and mole fractions z i , Equation (1) is solved to give n 2 through n N P . Then, the compositions of phases are obtained using i = 1, 2, · · · , N C j = 1, 2, · · · , N P (2) In this study, we introduce to simply Equation (1) into a set of fractional equations The values of ξ j i vary between −1 and ∞. Although Equation (4) appears simple, formulating a robust numerical solution is rather complex, because of the presence of N C hyperplanes on which Equation (4) becomes singular. Precarious application of the Newton method, without proper consideration of the hyperplanes, can lead to non-convergence or convergence toward non-physical solutions. Leibovici and Neoschil [26] recognized that the physical constraint on the mole fractions of the phases, 0 ≤ n j ≤ 1 and N ∑ j=2 n j < 1, defines a hypertetrahedron in the n 2 , . . . , n N P space, and this hypertetrahedron cannot be dissected by any of the hyperplanes. This property ensures that F j in Equation (4) are continuously differentiable in the hypervolume that encloses the hypertetrahedron of the physically admissible solutions and the immediately adjacent negative flash region. By starting the initial guess within the hypervolume defined by Equation (6) and relaxing the Newton steps such that they do not cross the singular hyperplanes, Leibovici and Neoschil [26] were able to ensure that their Newton method does not generate a solution outside of Equation (6). However, the condition of F j being continuously differentiable does not guarantee that the Jacobian that is needed by the Newton method is always well conditioned. When Newton iteration approaches the boundaries of Equation (6), the Jacobian can become nearly singular, causing the Newton method to fail [29]. The abovementioned limitation of the Newton method motivated us to seek a method that does not require the computation of Jacobian or any derivatives of F j . Our method began with the recognition that since F j must be monotonically decreasing in the direction of increasing n j . Because of this property, when n j is increased to approach a singular hyperplane while all the other variables are kept constant, F j → −∞ . When n j is decreased to approach a singular hyperplane while all the other variables are kept constant, F j → +∞ . Additionally, when n j → +∞ while all other variables are kept finite, F j decreases and approaches zero; when n j → −∞ while all other variables are kept finite, F j increases and approaches zero. Figure 1 shows, qualitatively for a three-phase system, the monotonicity of F 2 and F 3 to n 2 and n 3 . The monotonicity of F j to n j , but not the other variables motivated us to design a successive-substitution method to solve Equation (4), as follows. First, an initial guess is made within the hypertetrahedron of physically admissible solutions, n 2 0 = n 3 0 = · · · = n N P 0 = 1/N P . Here, the subscripts to n j denote the number of iterative substitutions. We begin the first iteration by holding the values of n 3 0 through n N P 0 constant and seek n 2 1 to satisfy F 2 n 2 1 , n 3 0 , · · · , n N P 0 = 0. Once n 2 1 is identified, we hold n 2 1 and n 4 0 through n N P 0 constant and seek n 3 1 to satisfy F 3 n 2 1 , n 3 1 , n 4 0 , · · · , n N P This process continues till F N P = 0 is solved and subscripts of all n j are updated to level "1", which finishes the first iteration. The second iteration begins by holding the values of n 3 1 through n N P 1 constant and seek n 2 2 to satisfy F 2 n 2 2 , n 3 1 , · · · , n N P 1 = 0, and ends when all the subscripts of n j are updated to level "2". Iteration would then begin at level "3", and continue till all n j are converged. With this successive-substitution strategy, solution of Equation (4) is reduced to successive solutions of single-variable equation F j n j = 0. Recognizing that F j n j is monotonic between its poles, we used the following method to solve F j n j = 0: • Evaluate F j n 2 m , · · · , n j−1 m , n j * m−1 , n j+1 m−1 , · · · , n N P m−1 . Here, m refers to the level of iteration and n j * in the square bracket is the variable that needs to be updated to level "m" such that F j n 2 m , · · · , n j−1 m , n j m , n j+1 m−1 , · · · , n N P m−1 = 0. • Determine the direction along which to vary n j * . If F j n 2 m , · · · , n j−1 m , n j * m−1 , n j+1 m−1 , · · · , n N P m−1 > 0, n j * should be increased. If, however, F j n 2 m , · · · , n j−1 m , n j * m−1 , n j+1 m−1 , · · · , n N P m−1 < 0, n j * should be decreased.

•
Check whether a solution exists along the direction of increasing (or decreasing) n j * . This is achieved by performing a line search to see whether increasing (or decreasing) n j * would lead to an intersection with a singular hyperplane. If n j * needs to be increased and increasing n j * generates intersections with singular hyperplanes, one solution must exist between the current n j * and the nearest intersection, because, as n j * approaches the intersection, F j n j * approaches −∞. If n j * needs to be increased, but there is no intersection with singular hyperplanes along the direction of increasing n j * , there is no solution to F j n 2 m , · · · , n j−1 m , n j m , n j+1 m−1 , · · · , n N P m−1 = 0 and the calculation stops. Similar criteria apply when n j * needs to be decreased.

•
When solution exists, vary n j * to find n j * m . Because there is one and only one solution between n j * m−1 and the nearest intersection with singular hyperplanes, a bi-section method is used. Testing begins at the midpoint between n j * m−1 and the nearest intersection. The interval that contains the solution is continuously halved until F j evaluated at the midpoint of the interval becomes zero within a prescribed precision.
Successive substitution and bi-section methods that are presented above are very robust. As the initial guess is within the hypervolume defined by Equation (6) and solution search is bounded by the singular hyperplanes, this algorithm will never generate a solution outside the hypervolume. A flow chart for this algorithm is provided in Appendix B.
At the end of this section, we present an example that involves gas, oil, and water phases. The system that is considered here only consists of three components: methane, n-butane, and water. The pressure is 13.8 MPa and the temperature is 93.3 • C. Methane, n-butane, and moisture are present in the gas phase (g = 1); methane, n-butane, and dissolved water are present in the oil phase (o = 2); and, dissolved methane, dissolved n-butane, and water are present in the water phase (w = 3). The equilibrium ratios needed by the calculation were obtained from solubility data of pure hydrocarbons in pure water [33], published equilibrium ratios for hydrocarbons in a retrograde gas [34], moisture content in natural gas [35], and water solubility data in hydrocarbons [36]. At the pressure and temperature of interest, the mole fractions of methane and n-butane in the water phase are: Assuming that the equilibrium ratios for methane and n-butane, between gas and oil phases, are not affected by whether these phases contain water or not, we obtained K go C 1 = 2.181, K go nC 4 = 0.350. Data on moisture content in natural gas suggests that at the pressure and temperature of interest, the mole fraction of water in the gas phase should be x Correlation in Hibbard and Schalla [36] suggests that mole fraction of water in the oil phase can be set to x o H 2 O = 0.004. The calculated compositions of the three phases and the equilibrium ratios are listed in Table 1, together with their corresponding ξ j i . This example was solved using a mixture composition of z C 1 = 0.6, z nC 4 = 0.35 and z H 2 O = 0.05. The initial guess was ( n g , n o , n w ) = (1/3, 1/3, 1/3). The case converged within five iterations with convergence criterion [∆ n o , ∆ n w ] ≤ 10 −7 . Figure 2 shows the path of convergence from the initial guess. Solution for n g , n o , and n w is (0.6725, 0.2981, 0.0294).
We note that this SS-RR method can be used to obtain, for three-phase equilibrium calculations, the lines on which F 2 = 0 and F 3 = 0. These lines are of special interest, because their intersection is the needed solution. Such lines have been used by Haugen et al. [37] and Li and Firoozabadi [24] in order to explain their area-based bi-section method to solve F 2 and F 3 . The dotted lines in Figure 2 mark the locations where F 2 = 0 and F 3 = 0, respectively. They were constructed from the paths of convergence of 1780 SS-RR calculations, each with a different initial guess within the area that is defined by Equation (6). These calculations all converged to the same solution, which indicates that initial guess is not important as long as it is bounded by Equation (6).  Figure 2 shows the path of convergence from the initial guess. Solution for g n  , o n  , and w n  is (0.6725, 0.2981, 0.0294).  Table 1. Circle is the initial guess and cross is the converged solution. The smaller shaded triangle is the domain of the physically permissible solutions. The larger shaded triangle is the domain defined by Equation (6)  We note that this SS-RR method can be used to obtain, for three-phase equilibrium calculations, the lines on which Firoozabadi [24] in order to explain their area-based bi-section method to solve 2 F and 3 F . The dotted lines in Figure 2 mark the locations where 2 0 F = and 3 0 F = , respectively. They were constructed from the paths of convergence of 1780 SS-RR calculations, each with a different initial guess within the area that is defined by Equation (6). These calculations all converged to the same solution, which indicates that initial guess is not important as long as it is bounded by Equation (6).  Table 1. Circle is the initial guess and cross is the converged solution.
The smaller shaded triangle is the domain of the physically permissible solutions. The larger shaded triangle is the domain defined by Equation (6) including the negative flash region. The dashed lines are the singular lines. The dotted lines inside the larger shaded triangle correspond to F 2 = 0 and F 3 = 0, the intersection of which is the converged solution.

Hybrid Newton-Successive Substitution Method
The convergence speed of the above SS-RR method is slow when compared to the NM-RR method of Leibovici and Neochil [26]. However, as SS-RR is robust and it does not require the Jacobian or any derivative of F j , it can be integrated into NM-RR to handle situations with poorly conditioned Jacobians.
In hybrid Newton-Successive substitution (NSS-RR) method, the Newton step is identical to that in Leibovici and Neochil [26]. Prior to each Newton step, however, the condition number of the Jacobian matrix, κ, is evaluated. A Newton iteration is carried out if κ is less than a prescribed value and a successive-substitution iteration is carried out if otherwise. A flow chart for NSS-RR is provided in Appendix B.
When starting a success-substitution iteration, the result from the previous Newton step can be directly used. When starting a Newton step after a successive-substitution iteration, however, care should be taken because at the end of successive substitution F N P = 0. Assume that n 2 m−1 , · · · , n N P m−1 is the result of the m − 1-th iteration and n 2 m , · · · , n N P m is the result of the m-th iteration, which is a successive-substitution. If the m + 1-th iteration is a Newton step, we recommend not to use n 2 m , · · · , n N P m to start the Newton step, but a combination of n 2 m−1 , · · · , n N P m−1 and n 2 m , · · · , n N P m . We specifically used to start the Newton step after successive-substitution. Geometrically, n i * is the average of the points along the path of successive-substitution from n 2 m , n 3 m−1 , · · · , n N P m−1 to n 2 m , n 3 m , · · · , n N P m . As each of these points makes one and only one F j zero, n j * calculated from Equation (8) is guaranteed to separate from surfaces F j = 0, thus making it a good choice in our opinion to start the Newton step.
In what follows, we present two examples for which NM-RR of Leibovici and Neoschil [26] did not converge to the correct answers, whereas both SS-RR and NSS-RR were successful. The first example is a fifteen-component, three-phase mixture. The second example is a twenty-component, five-phase mixture. In NM-RR and NSS-RR, the relaxation parameter that regulates the Newton step when it intersects with the singular hyperplanes was set to 0.5. The maximum condition number in NSS-RR was set to 10 10 . In all of the methods, the convergence criterion set on the difference between vectors n i m−1 and n i m was 10 −7 . Note that if the convergence criterion is tightened, then NM-RR can find the correct solutions for these examples. We are only using these examples to compare NM-RR and NSS-RR under the same convergence criterion. Table 2 gives the composition of the fifteen-component, three-phase mixture, the equilibrium ratios, and the equilibrium composition of the three phases. For this example, both NSS-RR and SS-RR converged to the correct root (−0.01686263294, −1.1254155641). NM-RR of Leibovici and Neoschil [26] however converged to a wrong root (−0.04078420653, −1.1004615900). As shown in Figure 3, many Newton steps were carried out very close to a singular line of the mixture. At the 21st step, a very large condition number of 5.0977 × 10 11 was encountered. The hybrid method carried out a successive-substitution step at this location and the subsequent Newton steps converged to the correct solution. The NM-RR method, on the other hand, lost accuracy at this location and it converged to a wrong solution. Figure 4 shows the variations in the condition number of the Jacobian during iterations. NM-RR finished in 29 iteration steps and NSS-RR finished in 28 iteration steps. The 21st step in NSS-RR is the only successive-substitution step carried out. Table 3 gives the composition of the twenty-component, five-phase mixture, the equilibrium ratios, and the equilibrium compositions of the five phases. For this example, both NSS-RR and SS-RR converged to the correct root (−0.00538660799, −0.00373696250, −0.00496311432, −0.00415370309). NM-RR of Leibovici and Neoschil [26], however, converged to a wrong root (−0.00287415017, −0.00392609623, −0.00798417906, −0.00350187286). Figure 5 shows that, at the 38th step, the condition number of the Jacobian reached 3.52 × 10 11 . Upon detecting this large condition number, NSS-RR performed a single successive-substitution step, which helped the following Newton steps to converge to the correct solution using a total of 54 iterations. NM-RR of Leibovici and Neoschil [26], however, lost accuracy after the 38th step and converged to a wrong solution.    As observed from the above examples, as SS-RR does not use any Jacobians nor derivatives, it is a useful method to switch to upon detection of a poorly conditioned Jacobian. In our hybrid NSS-RR, successive substitution is only activated on rare occasions, and hence it does not add significantly to the computational time. Checking the condition numbers of Jacobians, on the other hand, generated a non-negligible overhead to the algorithm. In average, the run time of NSS-RR is about 1.4 times of that of NM-RR.   Table 3. Composition and equilibrium ratios of the 20-component, 5-phase example.

Validations and A Three-Phase Equilibrium Calculation
After proving the robustness of the developed SS-RR and NSS-RR, in this section, we first compare our solutions of three-phase Rachford-Rice equations to those that are available in the literature. Then, we combine the developed solver with equation of state and Henry's law to solve a complete oil-gas-water three-phase equilibrium.

Validations of Three-Phase Solutions
Three cases of three-phase equilibria from Nichita et al. [38] are used here to validate our solutions of three-phase Rachford-Rice equations. The first case is a ternary mixture containing CO 2 , CH 4 , and normal-hexadecane (nC 16 ). The second case is a sour gas system with six components, the details of which can be found in Robinson et al. [39]. The third case is a quaternary mixture, the phase behavior of which was studied by Kohse and Heidemann [40]. Tables 4-6 list the compositions of the mixtures and the equilibrium compositions of the three phases. Table 4. Composition of CO 2 -CH 4 -nC 16 ternary mixture [38] and equilibrium compositions of the three phases (g: gas; l 1 and l 2 : liquids).  Table 5. Composition of the sour gas system mixture [38,39], and equilibrium compositions of the three phases (g: gas; l 1 and l 2 : liquids).  Table 6. Composition of the quaternary mixture [38,40] and equilibrium compositions of the three phases (g: gas; l 1 and l 2 : liquids). Equilibrium ratios that were obtained from these tables were passed to our multiphase Rachford-Rice equation solvers. We note here that all solvers, NM-RR, SS-RR, and hybrid NSS-RR, converged on these cases and gave identical results. In Table 7, we compare our results to those that were reported in Nichita et al. [38]. This comparison shows that results from our algorithm are in very good agreement with previously reported solutions. Table 7. Mole fractions of the three phases (g: gas; l 1 and l 2 : liquids) from this study and from Nichita et al. [38]. In a recent work by Okuno et al. [29], four three-phase examples were presented ( Table 1 in reference [29]), including flash and negative flash near critical points. We note that our solvers can achieve convergence on all four examples. Although a quantitative comparison is not possible because the authors did not report the final mole fractions, our converged solutions do visually agree with that are those presented in the figures of reference [29].

Three-Phase Equilibrium of An Oil-Gas-Water System
In this section, we combine multiphase Rachford-Rice solvers with equation of state (EOS) and Henry's law, and use the SS-PS approach to solve an oil-gas-water equilibrium. Water exists extensively in many hydrocarbon reservoirs. Water in reservoir may come from initial water traps, migration from nearby aquifers, or water injection for improved/enhanced oil recovery. The presence of water can affect the phase behavior of reservoir fluids, because some gases, in particular, CO 2 and H 2 S, are soluble in the aqueous phase. Three-phase equilibrium calculations are needed to correctly describe the phase behaviors of such systems [41,42].
In this calculation, fugacities of components in the vapor and liquid phases were modeled using Peng-Robinson EOS following standard procedures. Fugacities of components in the water-rich phase were modeled using Henry's law [43]. Henry's law for a component sparingly soluble in the aqueous phase is where H i is Henry's law constant of component i in the aqueous phase. Details of Henry's law can be found in Appendix A. The fugacity of water in the aqueous phase can be obtained from the fugacity of the solutes using the Gibbs-Duhem equation, as in Prausnitz [44].
Φ ws is the fugacity coefficient of pure water at the saturation vapor pressure p ws vp and v m is the molar volume of pure water. The fugacity coefficient of water Φ ws is calculated by the Chou equation, reported by Rowe and Chou [45], as follows: For a three-phase system with n c components, the following equations must be satisfied: Li and Nghiem [46] defined the following sets of equilibrium ratios that are consistent with our equilibrium ratios, as defined in Section 2. The three-phase Rachford-Rice equations for our vapor-liquid-water system, following Equation (1), are The equilibrium calculation procedure began by passing initial guesses of equilibrium ratios to the Rachford-Rice equations. n l and n w solved from the Rachford-Rice equations were used to compute compositions of the phases: x v i , x l i , and x w i , following Equations (4) and (5). These compositions were then used to determine the fugacity coefficients Φ v i and Φ l i and the fugacity f w i , from which the equilibrium ratios were updated. The equilibrium ratios were again passed to the Rachford-Rice equations, until the results converge.
We applied the method that is presented above to a six-component mixture, the composition of which is shown in Table 8. Detailed description of this case can be referred to Li and Nghiem [46]. We conducted equilibrium calculation at 10 MPa and 100 • C. Our calculated results are in good agreement with results in the literature, as shown in Table 9.  Table 9. Results for the case with three phases and six components in presence of water.

This Study
Li and Nghiem (1986) Two-Phase To illustrate the difference between a three-phase equilibrium calculation and two-phase equilibrium calculations that discount the influence of the aqueous phase, we conducted a two-phase equilibrium calculation where water is removed from the mixture. The mole fractions of the non-aqueous components are reconstructed based on the relative fractions of them in the original mixture. The results of this two-phase calculation are also shown in Table 9. Comparison shows that although the two-phase approach neglecting the water component generated phase compositions that are similar to those from the three-phase calculation, the ratio between liquid and gas mole fractions, however, is noticeably higher. At equilibrium, the aqueous phase contains 14% of H 2 S and 2% of CO 2 by mole, and is therefore highly corrosive. This corrosive nature can only be captured by a three-phase equilibrium calculation.

Conclusions
This paper presents a successive-substitution method to solve multiphase Rachford-Rice equations. In this algorithm, Rachford-Rice equations are solved sequentially, and the solution of each individual Rachford-Rice equation is achieved using the method of bi-section. The direction of bi-section is determined by the monotonicity of the equation, whereas the limit of bi-section is determined by the locations of the poles. These considerations ensure that the algorithm always reliably converges to either a physically admissible solution or a solution in the adjacent negative flash region.
In compositional reservoir simulations, the number of phase equilibrium calculations performed is very large. Accuracy and robustness without a significant penalty on efficiency are therefore very important. A hybrid method was developed in order to combine the advantages of this successive-substitution method and the Newton method. In this hybrid method, a successive substitution step replaces a Newton step when the local Jacobian becomes poorly conditioned. This hybrid method can effectively suppress the errors in the Newton method of Leibovici and Neochil [26] owing to poorly conditioned Jacobian. In terms of computational time, successive substitution steps add very little to the overall computational cost. The computational time of this hybrid method is 1.4 times of that of the Newton method of Leibovici and Neochil [26], and this difference is mainly due to the overhead needed to check the condition numbers.
We presented seven examples to show the characteristics and the accuracy of the algorithm. The first example shows the convergence behavior of the successive-substitution method. The second and the third examples present the advantages of the hybrid method over the Newton method of Leibovici and Neochil [26]. In the rest of the examples, we verified our algorithm using results from the literature. Lastly, we combined the developed multiphase Rachford-Rice solvers with equations of states for the oil and gas phases and Henry's law for the water phase and completed a three-phase flash calculation for an oil-gas-water system with sour gases. Our result again is in very good agreement with that previously reported. It shows specifically that a two-phase equilibrium calculation neglecting the dissolution of sour gases in the water phase is not a good approximation of a true three-phase flash calculation.
Author Contributions: R.G. performed the investigation, wrote the original manuscript and made the revisions. X.Y. and Z.L. are supervisors of this study.
Funding: R.G. and Z.L. thank the financial support from the National Science and Technology Major Project under Grant NO. 2017ZX05009005. R.G. specifically thanks the China Scholarship Council for sponsoring R.G.'s visit to the Colorado School of Mines.

Conflicts of Interest:
The authors declare no conflict of interest.

A(T)H(T)
functions used to compute the specific volume v C cohesive energy density of water p ws vp is the saturated vapor pressure of water at temperature T, v ws m is the molar volume of water at p ws vp and T, and h wo − h ws is the enthalpy departure of liquid water at p ws vp and T. The saturated vapor pressure of water p ws vp was calculated using the Buck equation [49] p ws vp = 0.61121 exp((18.678 − T 234. 5 )( T 257.14 + T )) (A7) T is in • C, and p ws vp is in psi. Equation (A7) is applicable to liquid water: T > 0 • C. When T < 0 • C, for ice we used the following equation The molar volume v m was estimated from the specific volume v , which, in turn, was from a correlation by Chou reported in Rowe and Chou [45].
v is the specific volume in cm 3 /g, T is the temperature in K, p is the absolute pressure in kgf/cm 2 , and w NaCl is the weight fraction of NaCl in solution. The temperature functions are,