Two―Phase Liquid Production Allocation in Multilayer Producing Wells Using Tempe

Received 14 March 2012； accepted 23 May 2012

Abstract

This work presents a methodology of allocating oil rate and associated water cut to each individual layer using temperature measurements and total surface production of oil and water.

This paper consists of two parts. In part one； an analytical forward model is proposed for wellbore temperature response under two―phase production in a multilayer geometry， using a nodal representation of the well. This model accounts for the formation geothermal gradient， steady―state oil―water flow in the wellbore， friction loss and Joule―Thomson effect in the wellbore， contrast in the thermal and physical properties of oil and water， wellbore heat losses due to unsteady heat conduction in the earth， and the mixing of the fluid streams of contrasting temperature.

The second part shows the application of the above solution by applying inversion techniques on temperature data coupled with forward model to allocate water and oil influx from producing layers. The inversion result is verified using a variety of commingled flow problems including a field case of a deviated well producing an oil―water mixture from two active completions. Inversion results seem to be robust within +/―15% provided the temperature contrast between the commingled layers is at least one order of magnitude greater than the resolution of the temperature measurements （e. g.， 4°C contrast for 0.1°C resolution）

Key words： Temperature measurements； Analytical forward model； Two―phase production； Multilayer； Inversion； Field case； Temperature contrast； Commingled layers

Daoud， A. （2012）. Two―Phase Liquid Production Allocation in Multilayer Producing Wells Using Temperature Measurements. Advances in Petroleum Exploration and Development， 3（2）， ―0. Available from： URL： http：///index.php/aped/article/view/j.aped.1925543820120302.281

DOI： http：///10.3968/j.aped.1925543820120302.281

In completions producing commingled streams through multiple formations， information about zonal rates and its change with time is important for routine well diagnosis. Back allocation of production rates helps in identifying zones of high water cut or gas―oil ratio and taking actions to improve well productivity.

Allocation of production rates to individual zones is conventionally done by running a production―logging suite， which through multiple sensing modules （velocity， phase holdup， pressure/temperature） localizes the non―productive and poorly productive intervals and measures relative contribution of different zones. As we move towards frontier fields （subsea， deepwater） post―drilling intervention for zonal allocation is generally ruled out based on operational economics. This is also the case for land/platform developments with highly deviated wells.

With the advent of distributed temperature sensors （DTS） （Al―Asimi et al.， 2002； Brown & Hartog， 2002） continuous measurement of temperature along the well can be obtained， and appropriate mathematical models coupled with data analysis and interpretation techniques opens a new direction for production allocation in complex or inaccessible wells.

Earlier attempts have used temperature logs to allocate production rates in producing wells （Curtis & Witterholt， 1973）. However these methods were applicable to single―phase production for long production times and with large separation between the zones. A recent work （Li et al.， 1999） does analyze how to obtain two―phase profiles in producing wells using both temperature logs and flowmeter data； however this method introduces certain approximations that limit its applicability to a wide range of production rates.

1. STATEMENT OF THE PROBLEM

The objective of this work is to develop a robust algorithm for two―phase oil―water production allocation in multilayer wells based on temperature measurements. This is achieved by， first， developing a forward model for temperature prediction in multilayer producing wells. Second， using inversion techniques to invert the temperature for rate allocation for a special case of two―layer producing wells. For the method to be robust an analytical model was developed from first principles， taking into account the essential physics of the problem. This model was compared against models commonly exist in the literature （Ramey， 1962； Sagar， Duty & Scmidt， 1991） for single layer wells and the inversion algorithm was tested on synthetic examples and applied to data from a real field example.

2. SOLUTION METHODOLOGY

An analytical solution is developed for the wellbore temperature response starting with the basic mass， momentum， and energy balance equations. First the solution for a single―layer problem under single―phase flow has been developed and compared with other models in the literature （Ramey， 1962； Sagar et al.， 1991） and with a numerical model （developed in a compositional simulator with a wellbore thermal option）. Then the solution is extended to two―layer production and finally to two―layer and two―phase production. Detailed derivation of the forward model can be obtained from other references （Daoud & Jalali， 2004； Rabie， Daoud， El―Tayeb & Abdel Dayem， 2010）. A brief description of the forward model development and the inversion algorithm is described in Appendix A.

3. APPLICATIONS

In this section， different synthetic cases are used to test the robustness of the forward and the inversion algorithm under two severe conditions： a）‘high’ and ‘low’ temperature contrast between the producing zones （13 F and 3 F）， b） imperfect knowledge of the input parameters required by the temperature forward model. Next， this approach will be demonstrated on a field case， where distributed temperature sensors are installed in a two―zone oil/water producing well.

3.1 Synthetic Cases

Two types of synthetic cases are used to test the robustness of the approach， the first type is for single―phase two―layer producing wells and the second is for two―phase two―layer producing wells.

3.1.1 Single―Phase Two―Layer Producing Wells

Nine cases have been generated to test the accuracy of rate allocation from two―zone producing wells， the first case is for perfect knowledge of the input parameter in the forward model and the rest are for imperfect knowledge. All the nine cases are tested under high temperature contrast （13 F） and low temperature contrast （3 F） between the producing zones. Table 1 shows the input data used to generate the true temperature for the nine cases， where the true oil rates from the lower and upper zone are 800 and 200 STB/D， respectively. In all the nine cases， we started our initial guess with zero rates， which is equivalent to using the geothermal temperature as initial guess. Table 2 and Figure 1 show the different cases used to test the accuracy of the rate allocation for the low and high temperature contrast between the producing zones. Case 9 shows a combined error of some of the parameters. Also， it is important to mention that only the points near the producing zones are used without imposing any constraints during the inversion. However for cases 8* and 9* the total rate has been introduced as constraint in the inversion algorithm.

From the results of these synthetic cases， it was found that：

（1） In all the cases that have error in model parameters （10%）， the rate allocation from the two producing zones shows error less than 10 % except for cases 8 and 9. Case 8， where errors in the temperature measurements have been generated randomly from Gaussian distribution of zero mean and 0.1 F standard deviation， while case 9 includes all the errors in the model parameters in one single case.

（2） Improvement in the rate allocation is obtained after including total rate as constraint in the objective function. Also， another observation is that case 9， which includes all the errors in the model parameters， shows better results compared to case 8. This might be due to the fact that some errors have compensating effects.

3.1.2 Two―Phase Two―Layer Producing Wells

Four different cases have been generated to test the accuracy of two―phase rate allocation from two―zone producing wells. The four cases studied are for perfect knowledge of input model parameters under low and high temperature contrast between the producing zones and with or without imposing total rate for each phase as constraint in the optimization algorithm. Parameters in Table 1 are used to generate the true temperature for the four cases （with additional information about the water phase with a specific gravity of 1， viscosity of 0.31 cp and specific heat capacity of 1 BTU/lbm―F）. The true rates for each phase from the lower and the upper producing zones， respectively， are as follows： qo1 = 300 STB/day， qo2 = 200 STB/day， qw1 = 350 STB/day， and qw2 = 300 STB/day. Similarly， we started our initial model with zero rates for the two producing zones. Figure 2 shows the results of the inversion for the four different cases.

From the results of the above study， it was observed that in spite of using perfect knowledge of the input model parameters， imposing the total rate of each producing phase as constraint does not improve the rate allocation if the temperature contrast between the zones is low. However， it does significantly improve the rate allocation when the temperature contrast between the layers is high. High temperature contrast may be qualified as temperature contrast that is one or two orders of magnitude greater than the resolution of the measurements.

3.2 Field Case

A test has been conducted on a deviated well producing from two active completions yielding about 4800 BOPD and 1500 BWPD. Fiber optic line has been installed in the completion to measure the temperature inside the wellbore and also information about the geothermal gradient has been provided. The temperature contrast between the two producing zones for this well is about 7 F.

From a quick inspection of the temperature data， it could be observed that not all the perforated interval from the upper zone is producing. Basically only a portion of the upper producing interval has a decline in temperature due to the mixing of the fluid of low temperature from the upper zone and high temperature from the lower zone. The remainder of the upper interval does not show this reduction in temperature.

Some forward model parameters are uncertain， so we add them as parameters in the inversion （within bounds） in addition to the zonal rates. Table 3 shows the constraints imposed for the model parameters and the rates in the inversion. The constraints used for the model parameters are taken from those commonly used in the literature （McCain， 1990； Popov， Pribnow， Sass， Williams， & Burkhardt， 1999）. We started our initial model of zero rates from each zone due to the lack of any prior information about the two phase zonal rates.

Figure 3 shows the comparison between the measured and the modeled temperature after the inversion. The rate allocations from lower and upper zone， respectively， were found to be as follows： qo1 = 2850 bbl/day， qo2 = 1986 bbl/day， qw1 = 1174 bbl/day， qw2 = 310 bbl/day. Figure 4 shows the absolute errors between the modeled and the measured temperature across the section of the well under study.

From Figures 3 and 4， it can be shown that the modeled temperature is in good agreement with the measured temperature with a maximum difference less than 1 F. As seen from Figure 4 the higher error is in between the two producing zones， which might be due to the presence of specific completion devices that are not accounted for in the model formulation.

During this test the productivity of each zone was tested separately （zones could be switched on and off with down hole valves） and it was found that the gross production from the lower zone is much higher than for the upper zone for both oil and water. This was also indicated by the inversion algorithm. Figure 5 shows the error in zonal rate allocation for each phase and for the gross liquid production compared to the zonal productivity tests at final iteration. It was found that the rate allocation error ranges from 7% to 36% at the final iteration taking into consideration that we started with zero zonal rates for each phase as initial guess.

（1）?A forward model of temperature has been developed from the basic principle of mass， momentum， and energy balance to describe the temperature profile in a single phase single layer producing well. The model is tested against established analytical models in the literature and a numerical model and the results show its good accuracy in the vicinity of the producing intervals.

（2）?The forward model has been extended to two zones two phase （oil―water） producing wells to be used as a forward tool in inverting for the zonal rate by knowing the temperature measurements along the wellbore.

（3）?An inverse modeling technique using the Generalized Reduced Gradient （GRG） optimization algorithm is used to invert for the zonal rate allocation by minimizing the difference between the measured temperature and the calculated ones from the forward model.

（4）?Several synthetic examples have been studied to test the accuracy of zonal rate allocation from temperature measurements in two―layer producing wells under different conditions and the results reveal the following：

a）?For single―phase liquid production with high temperature contrast between the producing zones， the zonal rates can be allocated with good accuracy without imposing the total rate as constraint during inversion.

b）?For single―phase liquid production with low temperature contrast between the producing zones， the zonal rates can be allocated with good accuracy if the total rate is added as constraint during inversion.

c）?For two―phase （oil―water） production with high temperature contrast between the producing zones， the zonal rates can be allocated with good accuracy only after imposing the total rate for each production phase as constraint during inversion.

d）?For two―phase （oil―water） production with low temperature contrast between the producing zones， the zonal rates are difficult to be allocated even if the total rate for each phase is added as constraint， the problem shows a high non―uniqueness and the optimization mainly depends on the starting guess for the zonal rates. So， prior information to select a good starting guess for the rates， or another constraint must be added to the problem to make it well posed.

e）?From error analysis study on the effect of different model parameters on the accuracy of zonal rate allocation， it was found that error in temperature measurements has the highest effect on the accuracy of zonal rate allocation. Error in temperature measurements （± 0.3 F） leads to more than 10% error in zonal rate allocation after imposing total rate as constraint in inversion for low temperature contrast between the producing zones. So， for accurate rate allocation， high temperature resolution is required.

（5）?Both the forward and the inversion are tested on a field data taken from a deviated well with 7 F temperature contrast between its two producing zones. The allocation algorithm attributes most of the oil and water production to the lower layer as confirmed by zonal tests， and reproduces the measured DTS profile within 1 F across the inter―layer and producing intervals.

ACKNOWLEDGMENTS

The author wishes to thank Schlumberger Well Completions and Productivity for sponsoring this study， Thang Bui for the Eclipse―300 numerical comparisons and Younes Jalali for his valuable advices during working in this project.

REFERENCES

[1] Al―Asimi， M.， Butler， G.， Brown， G.， Hartog， A.， Clancy， T.， Cosad， C.， Fitzgerald， J.， Navarro， J.， Gabb， A.， Ingham， J.， Kimminau， S.， Smith， J.， & Stephenson， K. （2002）. Advances in Well and Reservoir Surveillance. Schlumberger Oilfield Review， 14（4）， 14―35.

[2] Bird， R. B.， Stewart， W. E.， & Lightfoot， E. N. （2002）. Transport Phenomena （Second edition）. New York： John Wiley & Sons.

[3] Brown， G. A.， & Hartog， A. （2002）. Optical Fiber Sensors in Upstream Oil and Gas. Journal of Petroleum Technology， 54（11）， 63―65.

[4] Curtis， M. R.， & Witterholt， E. J. （1973）. Use of the Temperature Log for Determining Flow Rates in Producing Wells. Paper SPE 4637 presented at the 48th Annual fall meeting of the society of petroleum engineers of AIME， Las Vegas， Nevada， Sept. 30 ― Oct. 3， 1973.

[5] Dake， L. P. （2001）. Fundamentals of Reservoir Engineering， Amsterdam， The Netherlands： Elsevier.

[6] Daoud， A.， & Jalali， Y. （2004）. System and Method for Determining Flow Rates in a Well. Patent application published under the Patent Corporation Treaty （PCT） number PCT/IB2004/002639 and World Intellectual Property Organization number WO 2005/035943 A1.

[7] Hasan， A. R.， & Kabir， C. S. （1991）. Heat Transfer During Two―Phase Flow in Well Bores： Part I― Formation Temperature. Paper SPE presented at the 6th Annual Technical Conference and Exhibition， Dallas， Oct. 6―9， 1991.

[8] Lasdon， L. S.， Waren， A. D.， Jain， A.， & Ratner， M. （1978）. Design and Testing of a Generalized Reduced Gradient Code for Non Linear Programming. ACM Transaction on Mathematical Software， 4（1）， 34―50.

[9] Li， H.， Zhu， D.， Lake， L. W.， & Hill， A.D. （1999）. A New Method to Interpret Two―Phase Profiles From Temperature and Flowmeter Logs. Paper SPE 56793 presented at the 1999 SPE Annual Technical Conference and Exhibition， Houston， 3―6 October.

[10] McCain， W. D. （1990）. The Properties of Petroleum Fluids. Tulsa： PennWell Books.

[11] Nocedal， J.， & Wright， S. （1999）. Numerical Optimization. New York： Springer―Verlag.

[12] Popov， Y. A.， Pribnow， D.， Sass， J.， Williams， C.， & Burkhardt， H. （1999）. Characterization of Rock Thermal Conductivity by High―Resolution Optical Scanning. Geothermics， 28（2）， 253―276.

[13] Rabie， R.， Daoud， A.， El―Tayeb， E.， & Abdel Dayem， M. （2010）. Production Allocation in Multi―Layers Gas Producing Wells using Temperature Measurements. Paper SPE 139261 presented at SPE Latin American & Caribbean Petroleum Engineering Conference， Lima， Peru， 13 December.

[14] Ramey， H. J. Jr. （1962）. Well Bore Heat Transmission. Journal of Petroleum Technology， 14（4）， 427―435.

[15] Sagar， R.， Doty， D. R.， & Scmidt， Z. （1991）. Predicting Temperature Profiles in a Flowing Well. SPEPE， 6（4）， 441―448.

[16] Schlumberger. （2010）. Eclipse 100 and 300 Technical Description.

APPENDIX A

A.1 Forward Temperature Modeling for Single Phase Liquid， Single Layer Producing Wells

Figure A.1 shows the thermal nodal analysis used to develop a mathematical temperature model by determining the temperature at each node using mass， momentum and energy balance equations. Table A.1 shows the temperature nomenclature at each node presented in Figure A.1.

Figure A.1

Schematic Diagram Showing the Nodal Temperature Analysis in Tubing and in Formation for Single Layer Producing Well

Table A.1

Temperature Nomenclature at Each Node

Node Nomenclature

1 Bottom hole formation temperature calculated from the earth geothermal gradient， Teibh

2 Bottom hole flowing fluid temperature， Tfbh1

2` Flowing fluid temperature at the top of the producing zone， Tfbh2

3 Formation temperature at the well/earth interface in the non―producing zone

4 Formation temperature calculating from the earth geothermal gradient in the non producing zone， Tei

5 Flowing temperature， Tf， and it can be at any depth Z from the producing zone

I― Producing Zone （Node 1 and 2）

The temperature equation derived between node 1 and 2 based on energy balance equation is as follows：

（A.1）

The pressure drop can be calculated using Darcy’s equation by knowing reservoir rock and fluid properties and assuming a steady state flow （Dake， 2001）：

（A.2）

II― Node 2 2`

The producing zone between node 2 and 2` is divided into equal intervals， each interval producing equal rate. The number of divisions depends upon the number of temperature measurements in front of the producing zone between node 2 and 2`. By applying a macroscopic mass and energy balance （Bird et al.， 2002） due to the mixing of two streams， we can get the temperature at any position in the well bore in front of the producing zone using the following equation：

（A.3）

Where， i = 2， 3，……..， n （n is the number of temperature measurements in front of the producing zone）， taking into consideration that Tfbh1 is calculated from Eq. A.1 and A.2. Also， Tei at each interval should be corrected due to the pressure drop across the perforation using the same Eqs. A.1 and A.2 using Tei instead of Teibh.

As the fluids produced from each interval inside the producing zone have equal rate and equal specific heat capacity， Cp， so Eq. A.3 can be written as follows：

（A.4）

Accordingly， temperature at node 2` will be：

（A.5）

It should be noted that for no production， Tfbh1 and Tfbh2 are the geothermal temperatures.

III― Non Producing Zone （Node 4 ―3）

As the fluid produced， heat is transferred by convection inside the well bore and some of this heat is lost by conduction to the non―producing formation. Thus， inside the non―producing zone， the transport phenomenon is only the heat energy due to heat loss from the well bore to the non―producing zone. So the only balance equation required is the energy balance equation.

By applying the energy balance equation between node 3 and 4， the approximated solution （Hasan & Kabir， 1991） of temperature at node 3 can be obtained as follows：

for tD ≤ 1.5

for tD > 1.5

（A.6）

Where，

（A.7）

Th is the temperature at node 3 and Tei is the temperature at node 4 known from the geothermal gradient.

（A.8）

（A.9）

IV― Well Path （Node 2`― 5）

As the fluid proceeds from node 2` to 5， heat energy is transported by convection and also mass and momentum are transported due to the fluid flow. So， energy， mass， and momentum balance equations are applied between node 2` and 5. The derived temperature between node 2` and 5 for liquid phase flow is given as follows：

（A.10）

Where，

（A.11）

（A.12）

U is the overall heat transfer coefficient and can be calculated from Eq. A.13 as shown in （Bird et al.， 2002） under the following conditions given by （Sagar et al.， 1991）：

Thermal resistance of pipe and steel are negligible compared to the thermal resistance of the fluid in the tubing/casing annulus，

Radiation and convection coefficients are negligible and can be ignored

（A.13）

The friction loss coefficient， f， can be obtained as follows：

Where Re is the Reynolds number and can be obtained as follows：

（A.14）

Tei in Eq. A.10 is calculated by knowing the earth temperature at the bottom hole （Teibh）， which is a fixed， and the earth temperature gradient using the following equation：

（A.15）

Eq. A.10 can be converted to a dimensionless form using the following dimensionless parameters：

（A.16）

（A.17）

（A.18）

（A.19）

Eq. A.10 becomes：

（A.20）

The boundary condition used to solve the ordinary differential equation， Eq. A.20 is

（A.21）

This means that the temperature at the bottom hole is equal to Tfbh. The solution of Eq. A.20 using the boundary condition of Eq. A.21 gives the profile of the dimensionless temperature inside the producing well as function of the dimensionless depth， ZD

（A.22）

Eqs. A.16. and A.17 are used to convert the profile from the dimensionless domain to the real domain by knowing the fixed fluid temperature at the bottom hole of the well， Tfbh， and the depth of the well， L.

A.2 Model Comparison

The objective of this comparison is two folds； the first is to test the accuracy of our developed temperature model， Eq. A.22， against the numerical （Schlumberger， 2010） and the well known temperature modeling in the literature （Ramey， 1962； Sagar et al.， 1991）. The second is to test the accuracy of Eq. A.22 after neglecting the Joule―Thomson effect due to pressure drop across the perforation. The reason for this is to reduce the high uncertainty associated with the calculation of Joule―Thomson coefficient as its calculation required information about the pressure drop across the perforation or information about the flow regime， some reservoir properties like permeability， drainage radius， etc.， which are associated with high uncertainty.

The analytical temperature models used for comparison with the numerical are as follows：

（1）?“Ramey’s Model （Ramey， 1962）”.

（2）?“Sagar Model （Sagar et al.， 1991）”.

（3）?“Model 1”， which is the temperature model given by Eq. A.22.

（4）?“Model 2”， which is the temperature model given by Eq. A.22 after neglecting Joule―Thomson coefficient due to pressure drop across the perforation.

The numerical modeling is obtained by using ECLIPSE 300 （Schlumberger， 2010） under the thermal option to be able to get the temperature at specified node inside the producing well.

Table A.2 shows the well and the reservoir fluid and rock data used by the above four different models and the numerical model obtained from Eclipse 300 under the thermal option to do the comparison.

Table A.2

Input Data for the Numerical （Eclipse 300） and the Different Temperature Forward Models

Well Data Value Units

qo， oil Rate 1020 STB/D

qw， water rate 0 STB/D

qg， gas rate 306 MSCF/D

GLR 300 SCF/STB

Well Length， L 6792 ft

rti （inside tubing radius） 0.9075 in

rto （outside tubing radius） 1.1875 in

rci （inside casing radius） 2.506 in

rco （outside casing radius） 2.75 in

rwb （well bore radius） 3.75 in

Heat transfer coefficient between the well bore and the formation， U 156.55 BTU/D―ft2―F

Well roughness 0.001 ft

oil heat capacity， Cpo 0.485 BTU/lbm―F

oil API 58

Inclination to the horizontal， θ 90

Rock thermal Properties

Thermal rock conductivity （Ke） 33.6 BTU/D―ft―F

Heat capacity of rock 0.2115 BTU/lbm―F

Geothermal Gradient （GT） 0.0274 F/ft

Temperature of the earth at the bottom hole of the well （Teibh） 237.2 F

Reservoir Properties

Permeability （homogenous reservoir）， K 200 md

thickness of the producing layer， h 20 ft

drainage radius， re 0.62 Ft

oil formation volume factor， Bo 0.91 BBL/STB

Simulation Time， t 168 hr

Figure A.2 shows the comparison of the four models with respect to the simulation results； all the models show good agreement except for Sagar model. This is due to the correlation used by （Sagar et al.， 1991） where the example used might by beyond the database from which Sagar’s correlation was developed.

Figure A.2

Comparisons with the Numerical Simulation Model

Figure A.3 shows the point―by―point comparison with Eclipse results for all the models except Sagar’s model. This figure shows the plot of the error between each model and simulation results versus depth. It is seen from Figure A.3 that model 1 and 2 give less error compared to the other models very close to the producing intervals. At large distance from the producing intervals， models 1 and 2 show higher errors compared to the other models. The main reason that model 1 and 2 gave good results near the producing interval compared to the numerical simulation model is that gas holdup is still very small so the assumption of single phase liquid used in developing Eq. A.22 is still valid. Also， the assumption of constant friction loss is reasonable whenever the gas holdup is very small which is the case near the producing intervals. While at large distance from the producing intervals the gas holdup is increasing and the assumption of constant friction loss and single phase liquid production does not hold that is why model 1 and 2 are not giving good results at higher distance from the producing interval.

As seen from Figure A.3， that the difference between model 1 and 2 are small and could be neglected， so model 2 is selected in order to reduce the uncertainty of imperfect knowledge about the parameters used in calculating the Joule―Thomson coefficient due to pressure drop across the perforations.

Figure A.3

Error Comparison Between the Numerical Simulation and the Analytical Models

A.3 Temperature Forward Modeling for Single Phase Liquid， Two Layers Producing Wells

Figure A.4 shows a thermal nodal analysis sketch for two layers， single―phase liquid production wells. The only difference between the single layer and the two layers production is in the nodal analysis between nodes 5―5`， nodes 8―7， and nodes 5`―9， also there is a minor change between nodes 2`―5.

Figure A.4

Thermal Nodal Analysis for Two―Layer， Single Phase Producing Wells

I― Node 2`―5

The ordinary differential equation given across these nodes for the single layer production， Eq. A.20， is the same as that used for the two layers production however a more general boundary condition is used to solve it， which is：

（A.23）

This general boundary condition allow to handle the two or more layers as the temperature at node 5` should be corrected for the mixing of two streams and the change of the rate. This as if we are dealing with the well as consisting of different sections each has the same equation but different boundary condition depending upon the temperature of the previous section.

The solution of Eq. A.20 using the general boundary condition， Eq. A.23 is as follows：

（A.24）

Where， Tfdbh is the temperature of entry and zdbh is the depth measured from the bottom of the well at the entry level. Similarly， Eq. A.16 and Eq. A.17 are used to convert the dimensionless temperature profile obtained from Eq. A.24 to the real domain.

II―Node 5 ― 5`

The modeling between node 5 and 5` is very similar to that between node 2 and 2` for the single layer production， where both mass and energy balance are applied. Also the assumptions used between node 2 and 2` are the same as between node 5 and 5` except the last assumption where the heat capacity of the two streams are not the same and also the mixing rates are not equal.

Similarly， we divided the producing zone into equal intervals， each interval producing equal rate. By applying a macroscopic mass and energy balance （Bird et al.， 2002） due to the mixing of two streams from the upper producing zone and the total rate obtained from the lower zone， we can get the temperature at any interval inside the producing zone using the following：

（A.25）

Where， i = 1， 2，…..n （n is the number of divisions or the number of temperature measurements in front of the upper producing zone）， Tf（i）： is the temperature at each interval in front of the producing zone， Tf（0） is the well bore temperature at node 5， Cp2： is the specific heat capacity of the fluid in the upper producing zone， q1， q2 ： is the total production from the lower zone and upper zone respectively， Cp（i）： is the specific heat capacity at each interval in front of the producing zone and is calculating as a rate weighting average according to the following equation：

（A.26）

Where， i = 1， 2， …….，n and Cp（0） = Cp1 which is the specific heat capacity of the fluid produced from the lower zone， which is constant through the section between node 2` and 5.

At node 5`，

（A.27）

（A.28）

It should be noted that Tei in Eq. A. 25 could be obtained directly from the geothermal gradient after neglecting Joule―Thomson effect due to the pressure drop across the perforation from the upper zone.

III― Node 8 7

The same equation described between node 3 and 4 for single layer production， Eq. A.6， can be used between node 8 and 7， where the flow rate is the total rate from the two producing zones.

IV― Node 5` 9

Eq. A.24 can be used to describe the temperature profile between node 5` and 9 by using the total rate （q1 + q2） instead of q1. Also， Cp between node 5` and 9 is equal to Cp5` as calculated from Eq. A.28.

A.4 Temperature Forward Modeling for Two―Phase Liquid， Two Layers Producing Wells

Figure A.5 shows the nodal analysis sketch for two―phase two layers production， where the two phases are oil―water （liquid）.

The extension of the modeling to two―phase flow depends upon recalculating the parameters of the modeling for the two―phase flow between each node. In the non producing zone as there is no fluid flow， only heat energy flow， so the change from single phase to two phase flow will not affect the temperature modeling between nodes 3 and 4 and node 8 and 7. The following are the model parameter calculation between each node for two―phase liquid flow.

Figure A.5

Thermal Nodal Analysis for Two―Layer， Two Phase Producing Wells

I― Node 2―2`

The temperature modeling as shown from Eq. A.4 is independent on the type of phase flow， so there is no change in the temperature modeling for two phase flow between those nodes.

II― Node 2`― 5

Eq. A.24 is used to get the temperature distribution between these two nodes. The physical parameters required by Eq. A.24 are calculated as a rate weighted average as follows：

（A.29）

（A.30）

（A.31）

III― Node 5 5`

Eq. A.25 and A.26 are used to calculate the temperature between these nodes by substituting q1 by （qo1 + qw1 ） and q2 by （qo2 + qw2）， Cp2 is calculated from the following equation

（A.32）

IV― Node 5` 9

Same as between node 2`― 5， where Cp， ? ， ρ are calculated as follows：

（A.33）

（A.34）

（A.35）

The extension of the temperature forward modeling to multi―layers two―phase flow is typically follows the same steps for the extension from single to two layers.

A.5 Inversion Algorithm

The objective of this paper is to use temperature measurements along the wellbore and invert these temperatures to allocate the rate from the producing zones. So， this inversion required an optimization algorithm to find the independent parameters， which are the zonal rates that minimize the following objective function：

（A.36）

Subject to any constraints on m， where （m） is a vector of independent parameters mainly the zonal rates and the thermal and physical properties of the formation and the produced fluid.

Different optimization algorithms （Nocedal & Wright， 1999） can be used to invert the temperature measurements for rate allocation. Here we used a ready―made optimization algorithm， which is the “Generalized Reduced Gradient” （GRG） method. Details about the algorithm can be found elsewhere （Lasdon， Waren， Jain， & Ratner， 1978）. The inversion algorithm requires only the observed temperature and the temperature calculated from the forward model.