PRODUCTION LOGGING INVERSION BASED ON DAS/DTS

Abstract

A method of optimizing production of a hydrocarbon-containing reservoir by measuring low-frequency Distributed Acoustic Sensing (LFDAS) data in the well during a time period of constant flow and during a time period of no flow and during a time period of perturbation of flow and simultaneously measuring Distributed Temperature Sensing (DTS) data from the well during a time period of constant flow and during a time period of no flow and during a time period of perturbation of flow. An initial model of reservoir flow is provided using the LFDAS and DTS data; the LFDAS and DTS data inverted using Markov chain Monte Carlo method to provide an optimized reservoir model, and that optimized profile utilized to manage hydrocarbon production from the well and other asset wells.

Claims

1. A method of optimizing oil production, comprising: a) providing one or more fiber optic cables in a well in an oil reservoir; b) measuring low-frequency Distributed Acoustic Sensing (LFDAS) data in said well during a time period of constant flow and during a time period of no flow and during a period of thermal perturbations; c) measuring Distributed Temperature Sensing (DTS) data from said well during said time period of constant flow, during said time period of no flow, and during said period of thermal perturbations; d) providing an initial model of reservoir flow using said LFDAS data and said DTS data; e) inverting said LFDAS data and said DTS data using a Markov chain Monte Carlo method to provide a production profile for said well; and f) using said production profile to optimize said well and future wells in said reservoir and produce hydrocarbons therefrom.

2. The method of claim 1, wherein said thermal perturbations are created using a well heater.

3. The method of claim 1, wherein said thermal perturbations are created by opening and closing said well.

4. The method of claim 1, wherein said thermal perturbations are created using fluid injections into said well.

5. The method of claim 1, wherein said thermal perturbations are created using a fluid injection into said well, wherein said fluid is colder than said reservoir.

6. The method of claim 1, wherein said thermal perturbations are created using a fluid injection into said well, wherein said fluid is warmer than said reservoir.

7. The method of claim 1, using equations (1) or (2) or (3) plus equations (4), (5), and (6), or mathematical equivalents thereof: $\begin{array}{cc}\mathrm{vR}\frac{\mathrm{dT}}{\mathrm{dx}}=v\frac{\mathrm{dR}}{\mathrm{dx}}\left(T_p-T\right)+U\left(T_f-T\right)+\mathrm{vRG}\frac{\mathrm{dz}}{\mathrm{dx}},& \left(1\right)\end{array}$ where v is a maximum flow velocity at a heel of said well, R is a spatial distribution of flow velocity profile, normalized to a range between 0 to 1, T is a well fluid temperature during stable production, T.sub.p is a produced fluid temperature at each perforation location, T.sub.f is a surrounding formation temperature, U is an overall heat transfer coefficient between said surrounding formation and said well fluid, G is a pressure, volume temperature (PVT) coefficient, x is a measured depth of said well, z is a true depth of said well; $\begin{array}{cc}\frac{T_{i+1}-T_i}{\Delta x}=\frac{R_{i+1}-R_i}{\Delta x}\frac{T_{\mathrm{pi}}-T_i}{R_i}+\frac{U}{{\mathrm{vR}}_i}\left(T_{\mathrm{fi}}-T_i\right)+G\frac{z_{i+1}-z_i}{\Delta x}.& \left(2\right)\end{array}$ where i refers to an i-th perforation cluster; $\begin{array}{cc}\frac{\mathrm{dT}}{\mathrm{dx}}=\left(\frac{\mathrm{dR}}{\mathrm{dx}}+\frac{U}{v}\right)\frac{\delta T}{R}+G\frac{\mathrm{dz}}{\mathrm{dx}},& \left(3\right)\end{array}$ where δT is a difference between temperature measurements during shut-in and stable production periods; $\begin{array}{cc}R\left(x\right)=1-\left(1-R_{\mathrm{bot}}\right)\frac{\underset{i}{\overset{P_x<x}{.Math.}}{P}_i}{\underset{i=1}{\overset{N}{.Math.}}{P}_i},& \left(4\right)\end{array}$ wherein R(x) is normalized flow velocity profile, wherein P.sub.i is a normalized productivity of said i-th perforation cluster, P.sub.x is a perforation cluster location, and R.sub.bot is a normalized flow velocity at a toe end of a sensing section;
∈.sup.2=∥T−T′∥.sub.2+λ∥vR−v′∥.sub.2,  (5) wherein ∈.sup.2 is a penalty function, T′ is a measured flowing temperature in said well, and V′ is a measured flow velocity in said well;
∈.sup.2=∥T−T′∥.sub.2+λ∥R−R′∥.sub.2,  (6) where R′ is a normalized flow-velocity ratio measured by DAS.

8. The method of claim 1, using equations (2), (4), and (5) or (6).

9. The method of claim 1, wherein said one or more fiber optic cables are temporarily installed.

10. The method of claim 1, wherein said one or more fiber optic cables are permanently installed.

11. A method of optimizing oil production, said method comprising: a) providing one or more fiber optic cables operably coupled to one or more interrogators in a well in an oil reservoir; b) measuring low-frequency Distributed Acoustic Sensing (LFDAS) data in said well during a time period of constant flow and during a time period of no flow and during a time period of perturbation of flow, wherein said perturbation of flow is created by thermal perturbations; c) simultaneously measuring Distributed Temperature Sensing (DTS) data from said well during said time period of constant flow and during said time period of no flow and during said time period of perturbation of flow; d) using one or more of equations 1-6 or their mathematical equivalents: i) inverting U and G using a gradient-descent based optimization while holding Pi fixed; ii) inverting P.sub.i using a Markov chain Monte Carlo optimization; iii) repeating step i-ii) thousands of times to generate thousands of initial models; iv) randomly perturbing P.sub.i of an initial model to create a new model and retaining said new model as a final model if a penalty function ∈.sup.2 for said new model is smaller than that of said initial model, and otherwise abandoning said new model; v) repeating step iv for each of said initial models to generate thousands of final models; vi) statistically analyzing said thousands of final models to obtain production allocation results; and e) using said production allocation to optimize production of hydrocarbon from said well and future wells in said reservoir.

12. The method of claim 11, wherein LFDAS uses <0.1 Hz.

13. The method of claim 10, wherein said thermal perturbations are created with a well heater.

14. The method of claim 1, wherein said thermal perturbations are created by fluid injection into said well.

15. The method of claim 1, wherein said thermal perturbations are created by opening and closing said well.

16. A method of predicting hydrocarbon production allocation from a reservoir, said method comprising: a) providing one or more fiber optic cables operably coupled to one or more interrogators in a well in a reservoir, wherein one or more fiber optic cables are permanently installed in said well; b) measuring low-frequency Distributed Acoustic Sensing (LFDAS) data in said well during a time period of constant flow and during a time period of no flow, and during a time period of perturbations of flow; c) simultaneously measuring Distributed Temperature Sensing (DTS) data from said well during said time period of constant flow and during said time period of no flow and during said time period of perturbations of flow; d) using one or more of equations 1-6 or their mathematical equivalents: i) inverting U and G using a gradient-descent based optimization while holding Pi fixed; ii) inverting P.sub.i using a Markov chain Monte Carlo optimization; iii) repeating step i-ii) thousands of times to generate thousands of initial models; iv) randomly perturbing P.sub.i of an initial model to create a new model and retaining said new model as a final model if a penalty function ∈.sup.2 for said new model is smaller than that of said initial model, and otherwise abandoning said new model; v) repeating step iv for each of said initial models to generate thousands of final models; and vi) statistically analyzing said final models to obtain hydrocarbon production allocation results from said well.

15. The method of claim 16, wherein said perturbations of flow are created using a well heater.

16. The method of claim 16, wherein said perturbations of flow are created by opening and closing said well.

17. The method of claim 16, wherein said perturbations of flow are created by fluid injection into said well.

18. The method of claim 16, wherein said perturbations of flow are created by injecting a fluid into said well, said fluid being at a different temperature than said reservoir.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0092] FIG. 1 DTS measurements and well path for a hydraulically-fractured unconventional oil producer. The temperature is measured during shut-in and stable-production periods. Red crosses indicate the locations of perforation clusters. Dot-dashed line plots the true depth of the well.

[0093] FIG. 2 Flow velocity analysis using LFDAS signal. Left: DAS response during a well opening event. Right: Velocity estimation using semblance analysis.

[0094] FIG. 3 Two models with different production allocations but similar flowing temperature. Black dashed line plots the shut-in temperature, crosses show the perforation locations. Black solid line shows the measured flowing temperature, thick gray dashed line depicts the model prediction. Bar is the normalized productivity of perforation clusters.

[0095] FIG. 4 Results of the synthetic test with even productivity. Upper panel shows the temperature profiles and model predictions. Lower panel shows a statistical box plot of the final models. Dashed line shows the range of 5% and 95% of the models, while gray bars show 25% and 75%. Red square and green cross are mean and median of the models, respectively.

[0096] FIG. 5 Same as FIG. 4, except the true model is completion dependent.

[0097] FIG. 6 The averages of perforation cluster productivity per NCS. Left panel is for the even production model and right panel for the NCS-dependent production model. Dashed line defines 5% and 95% range of the final models, while the gray bar defines 25% and 75%.

[0098] FIG. 7 Inversion results from the real data example.

[0099] FIG. 8 Comparison between two different random sampling procedures to create initial models: a) the productivity of each perforation cluster is uniformly distributed, b) sampling procedure described in Section 5.1.

[0100] FIG. 9 Same as FIG. 4, except the input R.sub.bot is different from the true model and there is no constraint from the DAS velocity measurement (λ=0).

[0101] FIG. 10. Comparison between the flow velocity distribution R and DAS velocity measurements for the real data example in Section 4. Black lines are individual models and red cross shows normalized DAS velocity measurement: a) initial model distribution, b) final model distribution without velocity constrain (λ=0), c) final model distribution with velocity constrain (λ=1).

[0102] FIG. 11. An example of a sensing plan to acquire the DAS and DTS data for the production logging inversion.

[0103] FIG. 12. Same as FIG. 4 except a heater is placed at the end of the fiber to heat up the borehole fluid by 5° F.

[0104] FIG. 13. Schematic of well with two fiber optic cables therein connected to DAS and DTS interrogators/recorders.

DESCRIPTION OF EMBODIMENTS OF THE DISCLOSURE

[0105] Swan et al. (2017) developed a novel method that uses low-frequency DAS signal (LFDAS), which is sensitive to small temperature variations, to measure borehole flow velocities. This method was able to provide reliable flow velocity measurements for unconventional oil producers, but the results suffered from low spatial resolution and could not provide estimations at a spacing of perforation clusters.

[0106] Temperature measurements have long been used for production logging purposes (e.g., Ramey, 1962; Curtis, 1973). However, it is only applicable for conventional vertical wells, where geothermal gradients are substantial along the well path, and is not reliable for high-deviated horizontal wells due to the solution non-uniqueness (Ouyang, 2006).

[0107] By jointly inverting for DAS and DTS measurements, we can better constrain the production allocation with higher spatial resolution and associated uncertainty analysis. We have done this with Markov Chain Monte Carlo based inversion methods.

1 Statistical Procedures

1.1. MCMC-Based Methods

[0108] Markov chain Monte Carlo (MCMC) methods comprise a class of procedures used in statistics for sampling from a probability distribution. By constructing a Markov chain that has the desired distribution as its equilibrium distribution, one can obtain a sample of the desired distribution by observing the chain after a number of steps. The more steps there are, the more closely the distribution of the sample matches the actual desired distribution.

[0109] When an MCMC method is used for approximating a multi-dimensional integral, an ensemble of “walkers” move around randomly. At each point where a walker steps, the integrand value at that point is counted towards the integral. The walker then may make a number of tentative steps around the area, looking for a place with a reasonably high contribution to the integral to move into next.

[0110] Random-walk Monte Carlo methods make up a large subclass of Markov chain Monte Carlo methods. Random walk Monte Carlo methods are a kind of random simulation of Monte Carlo method. However, whereas the random samples of the integrand used in a conventional Monte Carlo integration are statistically independent, those used in Markov chain Monte Carlo methods are correlated. A Markov chain is constructed in such a way as to have the integrand as its equilibrium distribution.

[0111] Examples of random walk Monte Carlo methods include the following:

[0112] Metropolis—Hastings procedure: This method generates a random walk using a proposal density and a method for rejecting some of the proposed moves. This is sometimes called Metropolis-coupled Markov chain Monte Carlo (MCMCMC).

[0113] Gibbs sampling: This method requires all the conditional distributions of the target distribution to be sampled exactly. When drawing from the full-conditional distributions is not straightforward other samplers-within-Gibbs are used. Gibbs sampling is popular partly because it does not require any ‘tuning’.

[0114] Slice sampling: This method depends on the principle that one can sample from a distribution by sampling uniformly from the region under the plot of its density function. It alternates uniform sampling in the vertical direction with uniform sampling from the horizontal ‘slice’ defined by the current vertical position.

[0115] Multiple-try Metropolis: This method is a variation of the Metropolis—Hastings procedure that allows multiple trials at each point. By making it possible to take larger steps at each iteration, it helps address the curse of dimensionality.

[0116] Reversible-jump: This method is a variant of the Metropolis—Hastings procedure that allows proposals that change the dimensionality of the space. Markov chain Monte Carlo methods that change dimensionality have long been used in statistical physics applications, where for some problems a distribution that is a grand canonical ensemble is used (e.g., when the number of molecules in a box is variable). But the reversible-jump variant is useful when doing Markov chain Monte Carlo or Gibbs sampling over nonparametric Bayesian models such as those involving the Dirichlet process or Chinese restaurant process, where the number of mixing components/clusters/etc. is automatically inferred from the data.

[0117] Unlike most of the current Markov chain Monte Carlo methods that ignore the previous trials, using a new procedure called a Training-based Markov chain Monte Carlo or TBMCMC, the TBMCMC is able to use the previous steps and generate the next candidate. This training-based procedure is able to speed-up the MCMC procedure by an order of magnitude.

[0118] Interacting Markov chain Monte Carlo methodologies are a class of mean field particle methods for obtaining random samples from a sequence of probability distributions with an increasing level of sampling complexity. These probabilistic models include path space state models with increasing time horizon, posterior distributions w.r.t. sequence of partial observations, increasing constraint level sets for conditional distributions, decreasing temperature schedules associated with some Boltzmann-Gibbs distributions, and many others.

[0119] In principle, any Markov chain Monte Carlo sampler can be turned into an interacting Markov chain Monte Carlo sampler. These interacting Markov chain Monte Carlo samplers can be interpreted as a way to run in parallel a sequence of Markov chain Monte Carlo samplers. For instance, interacting simulated annealing procedures are based on independent Metropolis-Hastings moves interacting sequentially with a selection-resampling type mechanism. In contrast to traditional Markov chain Monte Carlo methods, the precision parameter of this class of interacting Markov chain Monte Carlo samplers is only related to the number of interacting Markov chain Monte Carlo samplers. These advanced particle methodologies belong to the class of Feynman-Kac particle models, also called Sequential Monte Carlo or particle filter methods in Bayesian inference and signal processing communities. Interacting Markov chain Monte Carlo methods can also be interpreted as a mutation-selection genetic particle procedure with Markov chain Monte Carlo mutations.

[0120] The advantage of low-discrepancy sequences in lieu of random numbers for simple independent Monte Carlo sampling is well known. This procedure, known as Quasi-Monte Carlo method (QMC), yields an integration error that decays at a superior rate to that obtained by IID sampling, by the Koksma-Hlawka inequality. Empirically it allows the reduction of both estimation error and convergence time by an order of magnitude.

[0121] More sophisticated methods use various ways of reducing the correlation between successive samples. These procedures may be harder to implement, but they usually exhibit faster convergence (i.e. fewer steps for an accurate result).

[0122] In addition to the above described MCMC based methods, new MCMC methods may be developed and used in the methods herein.

1.2 BFGS-Based Procedures

[0123] In numerical optimization, the Broyden—Fletcher—Goldfarb—Shanno or BFGS procedure is an iterative method for solving unconstrained nonlinear optimization problems. It belongs to quasi-Newton methods, a class of hill-climbing optimization techniques that seek a stationary point of a (preferably twice continuously differentiable) function. For such problems, a necessary condition for optimality is that the gradient be zero. Newton's method and the BFGS methods are not guaranteed to converge unless the function has a quadratic Taylor expansion near an optimum. However, BFGS has proven to have good performance even for non-smooth optimizations.

[0124] In quasi-Newton methods, the Hessian matrix of second derivatives doesn't need to be evaluated directly. Instead, the Hessian matrix is approximated using updates specified by gradient evaluations (or approximate gradient evaluations). Quasi-Newton methods are generalizations of the secant method to find the root of the first derivative for multidimensional problems. In multi-dimensional problems, the secant equation does not specify a unique solution, and quasi-Newton methods differ in how they constrain the solution. The BFGS method is one of the most popular members of this class. Also in common use is L-BFGS, which is a limited-memory version of BFGS that is particularly suited to problems with very large numbers of variables (e.g., >1000). The BFGS-B variant handles simple box constraints.

[0125] The above described BFGS-based and any newly developed variants thereof can be used in the methods described herein.

2. Logging Methodology

2.1 DAS Flow Velocity Estimation

[0126] DAS signal at very low-frequency band (<0.1 Hz) is very sensitive to small temperature perturbations, which can be used to track convectional thermal slugging during production (Swan 2017). For unconventional oil producers during stable production, the thermal slugging signals could be too small to be detected. As a result, extra well operations or tools have to be involved to create the required signal.

[0127] One method that can be applied on hydraulically-fractured unconventional producers is to repeatedly shut in and open the monitored wells, causing a pressure and/or temperature pulse to travel the well. After the hydraulic-fracturing operation, spatial variations of temperature in the formation near the borehole are created due to the uneven stimulation results (various injection volumes and fracture geometry at each perforation cluster). Thermal spatial gradients start to build up in the well bore during shut-in period through the conduction between the borehole fluid and the surrounding formation.

[0128] FIG. 1 shows an example of the measured temperature profiles during shut-in period (shut-in temperature) and stable-production period (flowing temperature) of an unconventional oil producer. It is obvious that perforation cluster locations are associated with the local minimums of the shut-in temperature due to the injection of cold fluid. As the well is opened after a period of shut-in, the borehole temperatures are perturbed due to the fluid convection, which can be detected by LFDAS. The signal can be used to estimate flow velocities. The detailed procedure for flow-velocity measurement is described by Swan et al. (2017) and Krueger et al. (2017).

[0129] FIG. 2 shows an example of the LFDAS signal and the corresponding velocity analysis result.

[0130] Another way to use LFDAS to measure flow velocities is to place a borehole heater at the end of the fiber. The heater is able to constantly perturb borehole temperature during stable production, which creates the thermal slugging signals in LFDAS for the velocity analysis. This method is more ideal because it directly measures flow velocity during stable production, but does require extra tool implementation. The measured flow velocities are served as inputs for the later inversion.

2.2 Temperature Based Model

[0131] The wellbore is simulated by a 1-D model that satisfies the conservation of mass, momentum, and energy. Ouyang (2006) presented a complete equation sets for a three-phase (oil, water, and gas) example. In this study, we demonstrate the procedure using a simplified single-phase model, which assumes different phases are well mixed in the borehole. This assumption is reasonable for fast-producing oil wells, and the simulation model can be replaced by a more complex multi-phase model if necessary.

[0132] A simplified 1-D wellbore temperature model during stable production can be presented as:

[00001]$\begin{array}{cc}\mathrm{vR}\frac{\mathrm{dT}}{\mathrm{dx}}=v\frac{\mathrm{dR}}{\mathrm{dx}}\left(T_p-T\right)+U\left(T_f-T\right)+\mathrm{vRG}\frac{\mathrm{dz}}{\mathrm{dx}},& \left(1\right)\end{array}$

where v is the maximum flow velocity at the heel, R represents the spatial distribution of flow velocities, normalized to the range between 0 to 1. T is the borehole fluid temperature during stable production, T.sub.p is the produced fluid temperature at each perforation location, T.sub.f is the surrounding formation temperature. U is a conductive heat transfer coefficient between the surrounding formation and the borehole fluid, which is determined by formation properties, well completion, as well as the heat capacity of fluid and phase combination. G is the PVT coefficient that describes a fluid temperature change when fluids pressure drop occurs due to lifting. Also, x is the coordinate along wellbore (measure depth) and z is true vertical depth of the wellbore.

[0133] This equation can be solved using a finite-difference approximation:

[00002]$\begin{array}{cc}\frac{T_{i+1}-T_i}{\Delta x}=\frac{R_{i+1}-R_i}{\Delta x}\frac{T_{\mathrm{pi}}-T_i}{R_i}+\frac{U}{{\mathrm{vR}}_i}\left(T_{\mathrm{fi}}-T_i\right)+G\frac{z_{i+1}-z_i}{\Delta x}.& \left(2\right)\end{array}$

[0134] For hydraulically-fractured unconventional oil wells, T.sub.p and T.sub.f can be approximated by the borehole temperatures measured after an extended shut-in period. This approximation simplifies EQ 1 to EQ 3:

[00003]$\begin{array}{cc}\frac{\mathrm{dT}}{\mathrm{dx}}=\left(\frac{\mathrm{dR}}{\mathrm{dx}}+\frac{U}{v}\right)\frac{\delta T}{R}+G\frac{\mathrm{dz}}{\mathrm{dx}},& \left(3\right)\end{array}$

[0135] where T is the difference between the temperature measurements during shut-in and stable production periods. This equation provides important insights for the later uncertainty discussion in Section 3.1.

2.3 Inversion Procedure

[0136] A direct gradient-based inversion of EQ 1 or EQ 3 leads to highly non-unique solutions which are initial-model dependent. FIG. 3 shows an example that two models with very different production allocations produce almost identical flowing temperatures using EQ 3. Because of this, we adapted an MCMC method to solve the problem stochastically.

[0137] U and G are treated as unknowns, since they are critical parameters in the model and are not usually constrained by laboratory results. U and G are not allowed to change spatially, assuming that there are no spatial variations of formation property, well completion, and fluid phase composition in the section of interest. Perforation cluster productivity is defined as a normalized number between 0 and 1. The flow velocity ratio R (x) can be presented as:

[00004]$\begin{array}{cc}R\left(x\right)=1-\left(1-R_{\mathrm{bot}}\right)\frac{\underset{i}{\overset{P_x<x}{.Math.}}{P}_i}{\underset{i=1}{\overset{N}{.Math.}}{P}_i},& \left(4\right)\end{array}$

[0138] where P.sub.i is the normalized productivity of i-th perforation cluster, P.sub.x is perforation cluster location, R.sub.bot is the normalized flow velocity at the end of the sensing section, which can be constrained by LFDAS velocity results and assumed to be known. The goal of the inversion is to find a combination of Pi, U, and G that fits both the flowing temperatures measured by DTS and the flow velocities measured by DAS.

[0139] The penalty function ∈.sup.2 is hence defined as:

∈.sup.2=∥T−T′∥.sub.2+λ∥vR−v′∥.sub.2,  (5) [0140] where T′ is the measured flowing temperatures, and V′ is measured flow velocities.

[0141] If the flow velocities are measured during transient periods using the method described in Section 2.1, the total production rate may be different from that during stable production. As a result, the penalty function should be defined as:

∈.sup.2=∥T−T′∥.sub.2+λ∥R−R′∥.sub.2,  (6) [0142] where R′ is the normalized flow-velocity ratio measured by DAS. In this step we assume the production allocation is the same between the transient period and stable production period.

[0143] We apply an iterative two-step inversion procedure to estimate the unknowns. The parameters U and G are inverted by using Broyden—Fletcher—Goldfarb—Shanno (BFGS) procedure, which is a standard gradient-descent based optimization (Byrd 1995).

[0144] The P.sub.i are inverted using MCMC, which is a random-walk based inversion procedure. For the P.sub.i inversion, starting from the initial model, P.sub.i is randomly perturbed at each iteration step to create a new model. If the value of the penalty function for the new model is smaller than that of the current model, the current model is updated. Otherwise, the new model is dropped, and the current model is randomly perturbed again. This process usually repeats thousands of times.

[0145] Thousands of initial models are created randomly (with the procedure described in Section 5.1) and perturbed to obtain a large set of final models. The productivity P.sub.i in each of the final models are then normalized so that the average productivity in each model is 1. Then the final models are statistically analyzed to obtain the production allocation results and uncertainties associated therewith. These results can then be used in continued production of oil, and/or in further optimization of the well and continued production therefrom.

[0146] The entire inversion procedure is described as follows: [0147] 1. Generate an initial model with randomized P.sub.i. [0148] 2. Invert for U and G using a gradient-decent based method while fixing Pi. [0149] 3. Invert for Pi using MCMC while fixing U and G. [0150] 4. Repeat step 2-3 multiple times to obtain a final model. [0151] 5. Repeat step 1-4 thousands of times to obtain a collection of final models. [0152] 6. Evaluate final models using statistical analysis and use that data in well or production optimization.

3. Synthetic Test

[0153] Two synthetic tests were designed to verify the inversion procedure. In both tests, we used the shut-in temperature in FIG. 1 as T.sub.p and T.sub.s. We set the normalized flow velocity R.sub.bot at the end of the sensing section to be 0.52, U to be 1×10.sup.−4 s.sup.−1, and G to be 7×10.sup.−4° F./ft. The synthetic flowing temperature was calculated using EQ 3. A uniformly distributed random noise with a peak amplitude of 0.1° F. was added. The inversion procedure was applied on the synthetic data and the inverted models were compared with the true model.

3.1 Even Production

[0154] FIG. 4 shows the inversion results for a synthetic model where the productivities of all perforation clusters are the same. One thousand final models were calculated using the proposed inversion procedure. Most of the models fit the data equally well, with the standard deviations of the predicted flowing temperatures being very small (<0.05° F.). The mean and median values of the perforation cluster productivity of the final models were very close to the true model, with an average error of 7.0% for the means and 7.7% for the medians. As used herein, the “true model” is the model used to generate the synthetic data for the inversion test. It is the true answer of the inversion result.

[0155] The uncertainty of the results—shown by dashed line and gray bar in the box plot— systematically increases towards to heel. This is due to the heel-ward increase of the total flow rate in the borehole, where the relative contribution of individual perforation cluster gets smaller. The uncertainty also increases with a smaller δT, because the flowing temperature is less sensitive to the perforation cluster productivity where the produced fluid temperature is the same or similar to the borehole temperature.

[0156] The inverted U has a mean value of 1.1×10.sup.−4 s.sup.−1 with a standard deviation of 7.2×10.sup.−5 s.sup.−1. The inverted G has a mean value of 4.4×10.sup.−4° F./ft with a standard deviation of 3.2×10.sup.−4° F./ft. The estimation of U and G can be improved if more vertical section of the well beyond the heel-most perforation is included. However, the included vertical section has to share the same formation thermal property and well completions as the horizontal section.

3.2 Completion Dependent Productivity

[0157] In this test, we mimicked a situation that there are three completion designs with different number of clusters per stage (NCS) existing in the sensing section. The designs alternated at each stage with the NCS being 7, 5, 3, 7, 5, 3, and 7 from the heel to the toe. We also assumed that in the true model, the productivity depended on NCS, with clusters in 3 NCS stages being 20% more productive than that in 5 NCS stages, and 50% more productive than that in 7 NCS stages.

[0158] FIG. 5 shows the results of this test. The fitting of the flowing temperature is equally good compared with the previous test, but the prediction of the true model is less accurate, with an average error of 15.4% for the means and 19.8% for the medians. The inverted U in this case has a mean value of 1:3×10.sup.−4 s.sup.−1 with a standard deviation of 7.7×10's.sup.−1. The inverted G has a mean value of 3.2×10.sup.−4° F./ft with a standard deviation of 3.1×10.sup.−4° F./ft.

[0159] Although the error for individual perforation clusters was substantial, the productivity differences from various completion designs can be clearly distinguished. The mean values of the perforation clusters with the same NCS were calculated for each final model and are summarized in FIG. 6. The dependency of the productivity on the completion designs can clearly be identified after the averaging of the results.

4 Real Data Result

[0160] We then applied the inversion procedure to real data collected from an unconventional oil producer. The well path, flowing and shut-in temperature profiles are as shown in FIG. 1. The fibers used for the DAS and DTS measurements were delivered into the borehole by a carbon rod, which only reached one-third of the horizontal producing section. The ratio of the total production from the section that is beyond the sensing section (R.sub.bot) is assumed to be 52%, which is constrained by the DAS flow velocity measurement (FIG. 2). Only the well section with measured depth deeper than 11300 ft was included in the inversion, due to a sudden change of geothermal gradient, which related to the conductivity of the formation rocks, at the measured depth around 11200 ft (FIG. 1). The results of the inversion are shown in FIG. 7. Unlike the synthetic test results, the results from the real data show some of the clusters producing statistically more than other clusters by a number of times.

5.1 Random Sampling Procedure

[0161] We found that a uniformly-distributed random productivity for each perforation cluster cannot efficiently sample all possible production profiles (FIG. 8A), because it has little probability of creating long-wavelength signals (a large well section that has larger or smaller productivity). We thus developed a new sampling procedure to generate long-wavelength signals.

[0162] Ten random numbers between 1 and R.sub.bot were generated and sorted with a descent order. Together with 1 and R.sub.bot, the twelve numbers were assigned as the R value at evenly spaced points within the sensing section, and linearly interpolated for the values in between. P.sub.i was then calculated using EQ 4. This method creates initial models that contain long-wavelength signals. The later random perturbation in the MCMC inversion added short-wavelength signals to the final models. FIG. 8 shows a comparison of R profiles between the uniform random sampling and our sampling procedure for 1000 initial models. Our sampling procedure covers a much larger range of possible flow rate profiles than the random sampling procedure.

5.2 DAS Velocity Constraint

[0163] For horizontal wells with small spatial temperature gradients, inversion results only based on temperature measurements are highly non-unique and could be biased. For demonstration purposes, we modified the synthetic test in Section 3.1. While keeping all the inputs the same, we changed the presumed bottom rate R.sub.bot to 70%, instead of 52% that the true model had. We also set the weighting parameter λ to zero to eliminate the constraint from the DAS velocity measurements.

[0164] The results, which are shown in FIG. 9, are significantly biased from the true model. On the other hand, if the bottom rate R.sub.bot is set correctly, the inversion results based only on temperature measurements should be consistent with the flow velocities which are measured independently by LFDAS.

[0165] FIG. 10 shows the comparison between the flow velocity distribution R, from the inversion results where the velocity constraint is eliminated (λ=0), and the DAS measured velocities (normalized) for the real data example in Section 4. The difference between the initial models (FIG. 10A) and final models (FIG. 10B) demonstrates the constraint of the temperature measurements on the results. It is clear that even without constraint of the velocity measurements (except for the bottom rate), the results based only on temperature measurements are consistent with the flow velocities measured independently by LFDAS. However, the uncertainty of the results can be further decreased by adding the velocity constraint (FIG. 10C).

5.3 Sensing Plan

[0166] In order to efficiently acquire the temperature and the flow-velocity measurements, a proper sensing procedure should be planned before the data acquisition. FIG. 11 shows an example of the sensing plan. The sensing plan should contain three key components: a long stable production period to stabilize the production flow and measure flowing temperature, an extended shut-in period to measure the formation temperature, and a transient period with a series of opening and closing operations or fluid injections to create the thermal slugging signals for the LFDAS flow-velocity measurements. If a borehole heater is available to create the required thermal perturbations, the transient period in the sensing plan can be removed, as the flow velocities can be measured during stable production.

5.4 Utilize a Borehole Heater

[0167] If a borehole heater can be deployed during DAS recording to create and increase the strength of the temperature signal, the production logging results can be significantly improved. First, the flow velocities can be directly measured by LFDAS during stable production, which are more reliable than the transient-period measurements. Secondly, if the heater is placed near the end of the sensing section, and is powerful enough to raise the borehole fluid temperature at the end of the fiber, the temperature difference between the flowing temperature and shut-in temperature can be artificially increased to lower the uncertainties of the inversion results.

[0168] FIG. 12 shows the inversion results of the NCS dependent production model in Section 3.2, except here the borehole fluid temperature at the fiber end was raised 5° F. by a heater. The heated flowing temperature and the original flowing temperature are inverted simultaneously. The errors of the resulted productivity are 7.5% for both means and medians, which are much smaller than the case without heater (>15%). Thus, we have shown that adding a heater to raise the temperature by about 5° F. significantly improved the reliability of the results.

[0169] The following references are expressly incorporated by reference in their entirety for all purposes. [0170] Bukhamsin, A., et al. (2016). Cointerpretation of distributed acoustic and temperature sensing for improved smart well inflow profiling. In SPE Western Regional Meeting. Society of Petroleum Engineers. [0171] Byrd, R. H., et al., (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16 (5), 1190-1208. [0172] Curtis, M., et al. (1973). Use of the temperature log for determining flow rates in producing wells. In Fall Meeting of the Society of Petroleum Engineers of AIME. Society of Petroleum Engineers. [0173] Dakin, J., (1985). Distributed optical fibre Raman temperature sensor using a semiconductor light source and detector. Electronics letters, 21 (13), 569-570. [0174] Hill, A. D. (1990). Production logging: theoretical and interpretive elements. [0175] Ouyang, L.-B., (2006). Flow profiling by distributed temperature sensor (DTS) system-expectation and reality. SPE Production & Operations, 21 (02), 269-281. [0176] Paleja, R., et al. (2015). Velocity tracking for flow monitoring and production profiling using distributed acoustic sensing. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. [0177] Ramey Jr, H., et al. (1962). Wellbore heat transmission. Journal of Petroleum Technology, 14 (04), 427-435. [0178] Van der Horst, et al. (2014). Fiber optic sensing for improved wellbore production surveillance. In IPTC 2014: International Petroleum Technology Conference. [0179] Vu-Hoang, D., et al. (2004). A novel approach to production logging in multiphase horizontal wells. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. [0180] US20140358444 Method of Hydraulic Fracture Identification Using Temperature [0181] US20170260842 Low Frequency Distributed Acoustic Sensing [0182] US20170260846 Measuring downhole temperature by combining DAS/DTS data [0183] US20170260849 DAS method of estimating fluid distribution [0184] US20170260854 Hydraulic fracture monitoring by low-frequency DAS [0185] US20170342814 Production Low-frequency DAS SNR improvement [0186] US20180045040 Production logs from distributed acoustic sensors [0187] US20180016890 Hydraulic fracture analysis [0188] U.S. Pat. No. 9,347,310 Multiphase flowmeter for subsea applications.