A Bayesian approach to deconvolution in well test analysis
AA Bayesian approach to deconvolution in well test analy-sis
Themistoklis Botsas, Jonathan A. Cumming and Ian H. Jermyn
Department of Mathematical Sciences, Durham University, Durham, UK.
E-mail: [email protected]
Summary .In petroleum well test analysis, deconvolution is used to obtain information about thereservoir system. This information is contained in the response function, which can beestimated by solving an inverse problem in the pressure and flow rate measurements. OurBayesian approach to this problem is based upon a parametric physical model of reservoirbehaviour, derived from the solution for fluid flow in a general class of reservoirs. Thispermits joint parametric Bayesian inference for both the reservoir parameters and the truepressure and rate values, which is essential due to the typical levels of observation error.Using a set of flexible priors for the reservoir parameters to restrict the solution space tophysical behaviours, samples from the posterior are generated using MCMC. Summariesand visualisations of the reservoir parameters’ posterior, response, and true pressure andrate values can be produced, interpreted, and model selection can be performed. Themethod is validated through a synthetic application, and applied to a field data set. Theresults are comparable to the state of the art solution, but through our method we gainaccess to system parameters, we can incorporate prior knowledge that excludes non-physical results, and we can quantify parameter uncertainty.
Keywords : Bayesian modelling, deconvolution, well test analysis, MCMC
1. Introduction
Within petroleum engineering, well test analysis encompasses a set of methodologies forthe planning, conducting, and interpretation of the results of controlled experiments on a r X i v : . [ s t a t . A P ] D ec T. Botsas, I. H. Jermyn and J. A. Cumming a well in a petroleum reservoir. The well test experiment controls the flow rate of fluidfrom the reservoir and observes the pressure at the bottom of the well (the ‘bottomholepressure’), producing a pair of time series. The goal of well test analysis is then to identifyand interpret the relationship between pressure and rate, and use this to make inferencesabout reservoir parameters, properties, geology, and geometry (Bourdet, 2002).A variety of technical tools have been developed in order to analyse the data fromwell tests (Gringarten et al., 2008), such as straight line analyses (Horner et al., 1951),type curve matching (Gringarten et al., 1979), and derivative curve analysis (Bourdetet al., 1983). However, these standard methods are typically restricted to the analysisof periods of constant-rate data, due to the corresponding simplification of the pressure-rate relationship under such conditions. Hence, variable rate data has to be treated asa collection of smaller independent data sets, resulting in inconsistencies in the analysisarising from the obvious dependency induced by the cumulative effect of production overtime.Deconvolution emerged as a major milestone in well test analysis (von Schroeter et al.,2004) due to its ability to deal with variable rate tests in their entirety, enabling theconsistent analysis of larger and richer data than was previously accessible. The abilityto deconvolve entire pressure and rate histories enables a coherent analysis of longer-term behaviour, providing greater insights into the flow behaviour of the reservoir andmore distant geological features.The key relationship for this inference is that, under appropriate regularity conditionsand ranges of reservoir conditions, the bottomhole pressure, ˜ p , and the flow rate, ˜ q ,written as functions of time t , are accurately related by:∆˜ p ( t ) = ˜ p − ˜ p ( t ) = (cid:90) ∞ g ( t − t (cid:48) ) ˜ q ( t (cid:48) ) dt (cid:48) = (cid:90) ∞ ˜ q ( t − t (cid:48) ) g ( t (cid:48) ) dt (cid:48) , (1)where ∆˜ p is the pressure drop from an initial equilibrium pressure ˜ p at t = 0, and usehas been made of the fact that g ( t ) = ˜ q ( t ) = 0 for t <
0, by causality and by definitionrespectively. The reservoir response function g is the object of primary interest, becauseit encapsulates all relevant well- and reservoir-specific information in a single function,and thus provides an effective summary and signature for the behaviour of a particularwell in response to production. The value g (∆ t ) describes the effect on the pressure at reprint t of the rate at t − ∆ t ; it is also the derivative of the pressure drop at time ∆ t inducedby a unit flow from t = 0. In principle, it can be derived as the coincidence limit ofthe Green function of the diffusion equation with ˜ q as a time-varying point source andwith diffusivity and mobility coefficients spatially-varying functions of reservoir and wellproperties.For visualisation, it is better to plot z ( τ ) = τ +ln g ( e τ ), where τ = ln t (Bourdet et al.,1983); Figure 1 shows an example. This greatly facilitates human interpretation of theresponse function, with features of the plot being directly associated with particular flowregimes and reservoir features. These features may be categorised by the times at whichthey appear, where time is effectively a surrogate for distance from the wellbore. • Early-time (a few seconds or minutes) is dominated by the wellbore and its im-mediate surroundings. The response usually begins with a unit slope straight line,representing the effect of fluid filling the uniform wellbore, typically followed by acharacteristic ‘bump’, which represents the impedance caused by localised damagedue to drilling, known as the ‘skin’. • Middle-time (a few hours or days) usually shows stabilization as the well drawsfluid from the relatively homogeneous region beyond the damaged skin, where a‘radial flow’ regime is established in which fluid moves towards the well from alldirections. • Late-time (weeks or months) indicates features of the reservoir far from the well-bore, such as impermeable boundaries or other heterogeneities.Figure 2 shows some common response shapes. Note that a specific shape may notuniquely identify a reservoir feature ( e.g. the fault response can also be achieved byappropriate changes in mobility) and that sequential features may mask each other.Of course, we do not observe entire pressure and rate functions: we only observevalues at a finite set of times. Furthermore, these observations do not give values of ˜ p or ˜ q directly, but rather values subject to measurement error and other uncertainties.von Schroeter et al. (2004) introduced the first deconvolution method that took theseissues into account using a form of penalized total least squares (TLS) regression, with T. Botsas, I. H. Jermyn and J. A. Cumming
Time P r ess u r e Time R a t e −3 −2 −1 0 1 2 − . − . . log Time Z + Fig. 1: The deconvolution of the bottomhole pressure, ˜ p (left), and the flow rate, ˜ q (middle) produce the reservoir response function g , whose transformation z (right) isplotted here.the response z represented as a piecewise linear function. The method has a numberof limitations. First, it requires, and is highly sensitive to, the specification of multiplehyperparameters; second, regularization is required to produce a smooth response func-tion; third, the flexibility of the linear spline representation means that the resultingresponse functions are not guaranteed to be physically possible; and fourth, it lacks acoherent approach to uncertainty analysis. To address these limitations, we propose a Bayesian approach. First, we replace thelinear spline representation of the response function by a parametric model based on anexplicit yet sufficiently general solution to the diffusion equation, thus avoiding the needfor regularization while restricting attention to physically sensible responses. Second, theaforementioned hyperparameters are treated as nuisance parameters and integrated out.Finally, the Bayesian approach inherently provides for a coherent uncertainty analysis.We first establish some notation. The well pressure, ˜ p , is observed at times t = { t i } ,idealised as point measurements, giving pressure data p = { p i } , for i = 1 , . . . , m . Wehave a single observation, p , for the initial reservoir pressure ˜ p . The rate is modelled asperiods of constant flow, ˜ q = { ˜ q j } , defined over known time intervals T = { [ T j , T j +1 ] } ,which are observed as q = { q j } , for j = 1 , . . . , N . We denote the piecewise constantfunction given by T and ˜ q , with zeroes elsewhere, as ˜ q .We are interested in the probability P ( g, ˜ q , ˜ p | p , q , p , κ ), where κ represents any reprint −2 −1 0 1 2 − . − . − . − . log Time Z (a) Infinite radial flow (unbounded reservoir). −2 −1 0 1 2 − . − . − . − . log Time Z (b) Single fault. −2 −1 0 1 2 − . − . − . − . log Time Z (c) Channel (two parallel faults). −2 −1 0 1 2 − . − . − . − . log Time Z (d) Closed reservoir (depletion). Fig. 2: Top view of indicative reservoir configurations (left subfigures) and resultingresponse functions (right subfigures).other prior knowledge e.g. T . Bearing in mind equation (1), we have P ( g, ˜ q , ˜ p | p , q , p , κ ) ∝ P ( p | ˜ p , g ∗ ˜ q, κ ) P ( q | ˜ q , κ ) × P ( p | ˜ p , κ ) P ( g | κ ) P (˜ q | κ ) P (˜ p | κ ) , (2)where ∗ indicates convolution, and we have used various independences to be describedlater. As stated, g will be given by a parametric model based on an explicit solution tothe diffusion equation, with parameters φ , so that the distributions on g will in fact beon φ . In the next section, we describe this model in detail.
2. The response model
The main challenge for a Bayesian approach to well test analysis is to identify an appro-priate prior distribution for the response function. The model should be flexible enoughto capture the shapes encountered in practice, restrictive enough to prohibit non-physicalresults, and computationally tractable. In addition, a model expressed in terms of welland reservoir properties would support more detailed interpretation than is possible witha generic representation such as a spline, permitting comparison to the results of otherwell test analyses; providing constraints on, and a physical interpretation for, the prior;and perhaps providing direct information about physical aspects of the system.
T. Botsas, I. H. Jermyn and J. A. Cumming
Here we use a variation of the ‘multi-region radial composite reservoir model’ (Acostaet al., 1994; Zhang et al., 2010) to satisfy these desiderata. First, through appropriatechoice of parameter values, the response function can be made to resemble the majorityof plausible response shapes. Second, since the model is derived as a solution to thediffusion equation representing the physical fluid flow problem, it is restricted to phys-ically sensible forms. Third, the finite-dimensional nature of the model makes it farmore computationally tractable than a description of an arbitrary reservoir. Finally, themodel itself is parameterized in a way that can be associated with the flow behaviour inthe reservoir (Bourdet, 2002).
R r
Region 1Region 2 r ω (a) Radial composite model. nn − i r ω R R R i R n − r (b) Multi-region radial composite model. Fig. 3: The radial composite reservoir models represent the reservoir as two (left) or more(right) concentric circular regions, each with different, constant mobility and diffusivityvalues, in the centre of which lies the wellbore.The multi-region radial composite model is based on the simple radial compositemodel with an infinite outer region (Satman et al., 1980). This model represents thereservoir as two concentric circular regions (the second region tends to infinity), at thecentre of which lies the wellbore, as shown in Figure 3a. The model has six parameters:three correspond to a well in a homogeneous reservoir (Bourdet et al., 1983), and affectearly-time behaviour and global scales; the remaining three characterize the transition reprint Parameter Name Scope Description Definition P Pressure match Global Pressure scale P = log (2 πh ( k/µ ) ) T Time match Global Time scale T = log (cid:16) πh ( k/µ ) C (cid:17) W Wellbore stor-age coefficient Early time Early time shape W = log (cid:16) Ce S πφc t hr w (cid:17) R Radius parame-ter Transition Dimensionlessdistance to tran-sition R = log (cid:16) r r w e − S (cid:17) M Mobility ratio Transition Ratio of mobilitybetween regions M = log (cid:16) ( k/µ ) ( k/µ ) (cid:17) η Diffusivity ratio Transition Ratio of diffusiv-ity between re-gions η = log (cid:16) ( k/µφc t ) ( k/µφc t ) (cid:17) Table 1: Parameters of the radial composite reservoir model.at the boundary of the inner region, and primarily affect the mid- to late-time behaviour.These parameters are sufficient to define the diffusion equation for the system, and henceto derive the response function. They are summarised in Table 1.The parameters can be written as functions of the physical properties of the wellboreand reservoir: the wellbore radius r w , the wellbore storage coefficient C , the wellbore skin S , the formation thickness h , the total system compressibility c t , the viscosity of the fluid µ , the permeability k and porosity φ of the medium, and the transition radius r . Therelationships between the derived model parameters and these fundamental parametersare given in the final column of Table 1. The parameters k , µ , φ , and c t are piecewiseconstant over the regions of the radial model. Changes in their values at the transitionboundaries represent reservoir features such as faults and changes in rock properties; dueto the resulting changes in the six derived parameters, these then manifest themselvesas identifiable features in the response function. However, the fundamental parametersare non-identifiable, and in the absence of strong prior information about their values,we will work here with the six derived parameters φ = ( P, T, W, R, M, η ). T. Botsas, I. H. Jermyn and J. A. Cumming
Figure 4 shows the effects on the features of the response function of changes in theparameters φ . The parameters T and P define time and pressure scales, and representthe responsiveness of the well and reservoir to production. The parameter W governs theimpedance to flow due to the skin effect surrounding the well. Parameter R correspondsto the distance from the wellbore to the inter-region transition, and so affects the time atwhich the impact of that transition is perceived. Finally, M and η describe the relativechanges in reservoir properties at the transition, resulting in a shift in the responsefunction stabilisation level, or a localised deviation from that level. In view of themultiplicative nature of the parameters T , P , W , M and η , logarithms are taken (base10 because this is the industry standard, in which values are most easily interpretable).The radius parameter R determines the time at which features appear in the responsefunction. In view of the similar spacing between response features when plotted on a logtime scale, and of the relationship between time and distance from wellbore, the radiusparameter is also treated on a log scale. −3 −2 −1 0 1 2 − . − . − . − . log Time z (a) T −3 −2 −1 0 1 2 − − − − log Time z −10123 (b) P −3 −2 −1 0 1 2 − . − . − . − . log Time z (c) W −3 −2 −1 0 1 2 − . − . − . − . log Time z (d) R −3 −2 −1 0 1 2 − . − . − . − . − . − . log Time z −1−0.500.51 (e) M −3 −2 −1 0 1 2 − . − . − . − . log Time z −1−0.500.51 (f) η Fig. 4: Response function sensitivity to changes in the parameters φ .The simple radial model can be extended to the multi-region case (Acosta et al.,1994; Zhang et al., 2010) by introducing additional concentric radial regions around the reprint model, resulting in a model of n regions and n − R i = log (cid:16) r i − r i − r w e − S (cid:17) , mobility ratio M i = log (cid:16) ( k/µ ) i ( k/µ ) i +1 (cid:17) , and diffusivityratio η i = log (cid:16) ( k/µφc t ) ( k/µφc t ) i +1 (cid:17) for each transition beyond the first. This gives a total of 3 n parameters: φ = ( P, T, W, { R i , M i , η i } ), with i = 1 , . . . , ( n − e.g. a shift of the stabilisation level, followed by a localiseddeviation from that level, requires at least two transitions, where, in the former there isa change in mobility, and in the latter, in diffusivity. However, it also raises the questionof how many regions to use, which is addressed in Section 5.
3. The Bayesian model
Having described the model for the response function g , we now describe the variouscomponents of equation (2). We first consider P ( p | ˜ p , g ∗ ˜ q , κ ) = P ( p | ˜ p , φ , ˜ q , κ ), which represents the observationalerror in the pressure data. We introduce a variance parameter σ p , and use a multivariateGaussian distribution for the observed pressures when σ p is known in addition to theother parameters: p | φ , ˜ y , σ p , κ ∼ N (cid:16) ˜ p ( φ , ˜ y ) , I σ p (cid:17) . (3)where I x denotes the diagonal matrix with x on the main diagonal, ˜ y = (˜ q , ˜ p ), ˜ p =˜ p − c ( g ( φ ) ∗ ˜ q ), and c samples any continuous function of t at the points t . The piecewiseconstant nature of ˜ q means that we can go further, and write c ( g ( φ ) ∗ ˜ q ) = C ( φ )˜ q ,where C ( φ ) is an ( m × N ) matrix; this makes the convolution computation particularlytractable. We note in passing that the time points t are not equally spaced, so that theobvious Fourier transform method for computing the convolution is not available.The magnitude of σ p is informed by the performance of the pressure gauges, which iswell-documented in the literature (Cumming et al., 2013a). For a typical well test, the T. Botsas, I. H. Jermyn and J. A. Cumming accuracy is estimated to be within ± σ p ∼ U(0 , σ p = 5.We now consider P ( q | ˜ q , κ ). Again, we introduce a variance parameter, and use amultivariate Gaussian model: q | ˜ q , σ q , κ ∼ N(˜ q , I σ q ) . (4)Expert judgement suggests an associated uncertainty in measured rates of up to 10%of their magnitude. We keep the value fixed to a constant σ q = 0 . q m , where q m =max j ( { q j } ).Turning to P ( p | ˜ p , κ ), we again introduce a variance parameter σ p and a Gaussianmodel: p | ˜ p , σ p , κ ∼ N(˜ p , σ p ) , (5)The standard deviation is usually fixed to σ p = 10, again informed by expert judgement(Cumming et al., 2013a). We now turn to the distributions P ( φ | κ ), P (˜ q | κ ), and P (˜ p | κ ). The independencebetween these quantities arises from their distinct physical meanings and the availableprior knowledge. φ Absent other information, we assume independence for each of the reservoir and systemparameters in φ , producing a distribution of the following form: P ( φ | κ ) = P ( P | κ ) P ( T | κ ) P ( W | κ ) n − (cid:89) i =1 P ( R i | κ ) P ( M i | κ ) P ( η i | κ )for a multi-region radial composite model with n regions. For each component of φ ,we approach prior specification from the perspective of eliminating non-physical andimplausible values as robustly as possible. Where possible, choices for priors are derivedfrom information in geological and geophysical studies of the corresponding parameters,supplemented by expert knowledge from petroleum engineers. Alternatively, when the reprint parameters of our model do not equate to quantities of particular geological interest, werely instead on a synthesis of expert judgement, inspection of results from the analysisof other well tests in the literature, and knowledge of the sensitivity of the responsefunction to the parameters as seen in Figure 4. Based on these considerations, we useGaussian priors T ∼ N(2 , . ) and P ∼ N(1 . , . ), thereby spanning around twoorders of magnitude in the (non-logarithmic) time and pressure scales; while for W , weuse W ∼ Ga(1 , . R i represent positive increments in distance from the previous boundary. Inorder to replicate our prior expectation of homogeneous behaviour in the early stages ofthe test, but also the assumption that there is reasonable distance between the differenttransitions, we use R i ∼ N(2 , i th region is 100 times the ‘effective wellbore radius’ r w e − S , while the variance means thatthe width varies between 0 . times the effective wellbore radius (up to threestandard deviations).For { M i } and { η i } , we adopt a prior mean value of 0, corresponding to no change be-tween regions, and a symmetric distribution, representing an absence of prior knowledgeabout whether these quantities will increase or decrease at a transition. To avoid degen-erate solutions in which the ratio values grow exceptionally large, we restrict ourselvesto a relatively modest variance, and use M i ∼ N(0 ,
1) and η i ∼ N(0 , variation in each direction. The individual true rate values are independent P (˜ q | κ ) = (cid:81) Nj =1 P (˜ q j | κ ). We use uniformimproper priors on R ≥ for the components: negative rates represent an injection of fluidinto the well, which does not occur in our data. Taking the sign of the rates into accountis important, since it might result to changes in the shape of the posterior responsefunction. We use the same prior for ˜ p . T. Botsas, I. H. Jermyn and J. A. Cumming
There is no closed-form expression for P ( φ , ˜ y , σ | p , q , p , κ ), where ˜ y = (˜ q , ˜ p ) as before,and σ = ( σ p , σ q , σ p ), due to the complex form of the response as a function of φ .Consequently, a numerical method such as MCMC is needed to extract informationfrom the posterior.We can simplify things somewhat, however, by factorizing the posterior: P ( φ , ˜ y , σ | x ) = P ( ˜ y | φ , σ , x ) P ( φ , σ | x ). We find a Gaussian conditional posterior for ˜ y : P ( ˜ y | φ , σ , x ) ∼ N( A − b, A − ) , (6)where A and b are given by A = mσ − p + σ − p − Tm I σ − p C ( φ ) − C ( φ ) T I σ − q m I σ − q + C ( φ ) T I σ − p C ( φ ) , and b = (cid:32) m (cid:88) i =1 p σ − p + p σ − p , q T I σ − q − p T I σ − p C ( φ ) (cid:33) T , where 1 a is the a -dimensional vector with all elements equal to 1. The marginal posteriorfor φ and σ is: P ( φ , σ | x ) ∝ | πA | exp (cid:26) − m ln ( σ p ) − N ln ( σ q ) − ln ( σ p ) − (cid:16) p T I σ − p p + q T I σ − q q + p σ − p (cid:17) + tr( b T A − b ) (cid:27) . Using this factorization, sampling from the posterior can be performed by sampling fromthe marginal posterior for φ and σ , and then sampling from the Gaussian for ˜ y . This hastwo advantages: first, that the dimensionality of the space for which MCMC is requiredis greatly reduced (in our examples, from around 200 to around 10); second, that thisdimensionality no longer depends on the length of the time series of measurements,which is critical for the scalability of the method. Even for modest sample sizes, theresult is a significant reduction in the time needed to generate a given number of samplesof ( φ , σ , ˜ y ); in the case that one is interested only in φ , the reduction is even moresignificant. reprint Sampling from P ( φ , σ | x ) requires a careful choice of MCMC algorithm. First, weexpect strong correlations between the posterior parameters: Gibbs-type samplers areunsuitable and all parameters need to be updated simultaneously. Second, the covariancestructure is not known a priori: the algorithm needs to be adaptive in order to discoverthe covariance structure as samples accumulate. Finally, we have reason to expect thatthe posterior will be multimodal, as different combinations of parameter values canyield similar response functions: the algorithm needs to be able to detect the presenceof multiple modes and jump between them.To tackle these challenges, we use the Differential Evolution Markov Chain withsnooker updater and fewer chains algorithm (DEzs), as described by ter Braak and Vrugt(2008). The premise is to construct an adaptive random walk Metropolis MCMC usingthe concept of differential evolution (DE) (Storn and Price, 1997). The MCMC analogueto DE is the Differential Evolution Markov Chain (DEMC) (ter Braak, 2006), where thecandidate solutions in DE correspond to N parallel chains. The proposal for each chainis derived from the remaining N −
4. Example: Synthetic channel
In this section, we validate the method by application to a synthetic reservoir with aknown response function with channel behaviour (Figure 5b). We generate a data set(Figure 5a) corresponding to a well test of length 320 hours, comprising three periodsof production at constant but different rate values, and a total of 272 pressure mea-surements generated by convolution of the rate with the known response. To represent T. Botsas, I. H. Jermyn and J. A. Cumming potential observational error, the pressure measurements include additive random noisewith known variance σ p = 25. Time0 50 100 150 200 250 300 P r e ss u r e R a t e (a) Rate and pressure measurements. −3 −2 −1 0 1 2 − . − . − . − . − . − . log Time Z (b) True response. Fig. 5: Synthetic channel data.Our primary focus with this analysis is verify that our methods are able to capture thereservoir response behaviour, even though a channel is not a radially composite reservoir,and to determine how precisely the parameter values can be obtained from the data givenonly weakly informative priors. Consequently, we treat the true rate values and initialpressure as known, and equal to their measured values, leaving φ as the parameter ofinterest, and we choose vaguer priors for φ than those recommended in Section 3.2: aGamma prior for W and broad uniform priors for the remaining parameters. The DEzsalgorithm was run with three parallel chains for 666,667 iterations. A thinning of 50 anda burn-in of 1,000 were applied, leaving 37,002 iterations to be used for inference.The resulting response function posterior samples are shown in Figure 6a, as lightblue curves; the MAP response is shown in dark blue and the true channel response isoverlaid in red. All samples show the characteristic limiting half-unit slope expectedof a channel reservoir, thereby validating the response shape. Figure 6b shows theresiduals between the posterior samples of reservoir pressure and the true values fromthe underlying synthetic model. Most of the samples lie within ± reprint (a) Posterior response functions. (b) Posterior pressure residuals.(c) Scatterplots. (d) Extended response functions. Fig. 6: Posterior marginal distributions, posterior response functions, pressure residuals,and extended response functions, for the posterior response parameters in the syntheticchannel model.pressures, which is well within usual tolerances for such analyses.Figure 6c shows summaries of the posterior MCMC samples of the response param-eters giving rise to the response curves in Figure 6a. The marginal posteriors for T , P , W , and R have substantially smaller variances than their priors, indicating that the dataidentifies these quantities quite precisely. In contrast, the marginal posteriors for M and η are generally close to their prior forms, indicating an inability to identify specific valuesfor these transition parameters. This picture is refined by the joint posterior marginals,shown as scatterplots in Figure 6c. First, we see that there is a strong positive correla-tion between P and W , and weaker positive correlations between both those parametersand T . These can be explained by the fact that small changes in the early-time shape T. Botsas, I. H. Jermyn and J. A. Cumming ( W ) can be replicated by small shifts in the response ( T and P ), e.g. larger P valuesshift the response curve downwards, while larger W values make the ‘bump’ steeper;these effects cancel each other over a narrow range.The late-time parameters ( R, M, η ) also show interesting behaviour. The (
M, η )plot shows that, unlike T , P , W , and R , these parameters cannot be well identifiedindividually, i.e. cannot be constrained to a quasi-0-dimensional subset, but they can beconstrained to a quasi-1-dimensional subset, and in fact, the relationship is essentiallyaffine. However, there are two separate regimes. For smaller R , where the bulk of theprobability lies, M and η can take on a large range of values (although themselves stilllinearly related); for larger R , M and η are essentially determined as a function of R .There is a clear physical explanation for these two regimes, best explained by referenceto Figure 6d. This shows the posterior response function samples from Figure 6a plottedover an extended range; the red curves correspond to the first regime, the blue to thesecond. In the first, the limiting slope is indeed present. Its time of onset (essentially R )is quite well specified, but provided the change in rock properties passes some threshold,the exact values of M and η do not matter. In the second, the limiting slope is anillusion produced by the window of observation; in reality, there is a jump in the responsefunction to a constant higher level. In this regime, the time of onset and the values or M and η are tightly coupled: the smaller the jump, the later must be the onset topush it out of the observation window. These results reveal that the standard industryinterpretation of the half-unit slope as exclusively indicative of a channel is too definitive:a range of other reservoir configurations are consistent with this behaviour.
5. Example: Field data from an oil reservoir
In this section, we analyse the results of a real oil field dataset of 2,273 pressure obser-vations and 22 production rates over approximately 150 hours, shown in Figure 7. Weconsider radial composite models with 1–4 transitions (2–5 regions). The true rates andinitial pressure are unknown, so that in addition to using the DEzs algorithm for φ and σ p , we sample from the Gaussian for ˜ q and ˜ p , as described in Section 3.3. We ran threeparallel chains in the DEzs algorithm for 1,666,667 iterations each, which, after thinning reprint Time0 50 100 150 P r e ss u r e R a t e Fig. 7: Field oil data rate and pressure measurements (a) Posterior response functions. (b) Posterior rates.
Fig. 8: Posterior response and rate samples for the 3-transition model.of 50 and burn-in of 10,000 were applied, resulted in 70,002 samples for inference.For the following discussion, we focus on the 3-transition model. The posterior re-sponse and rate samples are shown in Figure 8. The response shows a very early feature,at around τ = − .
5, suggesting that the wellbore is very close to an inhomogeneity. Therate samples in Figure 8b show significant deviation from the measured values, althoughthe general structure remains. This is consistent with the known scale of errors in thesemeasurements.The one- and two-dimensional marginals are shown in in Figure 9 as smooth densityestimates and scatterplots respectively. The marginals for T , P , and W tend to high T. Botsas, I. H. Jermyn and J. A. Cumming
Fig. 9: Scatterplots for the response parameters.values relative to their priors, but although their posterior variances are larger thanin the synthetic case, as one might expect, they are still small relative to the priors.The radii parameters are reasonably far apart, suggesting that there is no degeneracyarising from a trivial region. The ˜ p marginal suggests that the true initial pressureis larger than the corresponding measurement, but the difference, around 2psi, is notunreasonable. The variance σ p lies in the range [2 , . T ,which implies that the horizontal shift of the response has a more significant role incompensating for the effects of P and W . Specifically, larger values of T shift the curve reprint − − RD M (a) Uncertainty step plot for m . − − RD h (b) Uncertainty step plot for η . Fig. 10: Uncertainty step-plots for the 3 transition model applied to the oil field dataset.to the left, causing the skin ‘bump’ to appear earlier, which, in consequence, causes W to increase.There are also numerous late-time parameter correlations, some of which are rea-sonably interpretable. For example, we detect negative correlations between M andthe parameters M and R , indicating that larger first transitions correspond to smallersecond transitions. In terms of the response, this says that the earlier the effect of thesecond transition arrives, the larger the size of the transition needs to be to producesimilar results. Figure 10 shows an alternative representation of the posterior for the reservoir parame-ters, showing how mobility and diffusivity change with radial distance from the wellbore.Specifically, defining ρ i = log (cid:80) ij =1 e R j = log (cid:0) r i r w e − S (cid:1) and m i = (cid:80) i − j =1 M j = log (cid:0) ( k/µ ) ( k/µ ) i (cid:1) ,we plot the piecewise constant functions taking [ ρ i − , ρ i ) to m i and η i − respectively. Inthese plots, a solid line indicates the posterior mean, while the coloured bands depict theposterior sample range (dark blue), 99% credible posterior intervals (middle blue), and95% credible posterior intervals. In general terms, the step changes in these parameterscorrespond to the features in the response function after the early-time skin effect.The first change in parameter values consists of a small decrease in η and an increase in T. Botsas, I. H. Jermyn and J. A. Cumming (a) 1 transition model. (b) 2 transition model.(c) 3 transition model. (d) 4 transition model.
Fig. 11: Response uncertainty for 1, 2, 3 and 4 transitions respectively including theTLS result (as a single black curve). m : this produces the increase in the response curve at around τ = − m and the switch to positive values for η . Finally, the late-time unit slope observed inthe response function corresponds to the third, more substantial positive shift in m , andthe wide range of possible values of η . These changes are indicative of an impermeableboundary to the reservoir, and this is indeed the standard interpretation of a late-timeunit slope. In Figure 11 we plot the posterior response function samples for models with 1–4 transi-tions. Each plot includes a black curve showing the TLS result, the current state-of-the- reprint art. The 1-transition result looks somewhat different to and flatter than the other results,probably indicating an under-parameterization giving an over-smoothed result. The 2-transition result is closest to the TLS result, with slight deviations in early times, againprobably from the relative rigidity of the model. The 3- and 4-transition results are verysimilar to each other, suggesting that the 4-transition model is over-parameterized andthat the 3-transition model is an adequate description. The 3- and 4-transition resultsshow a qualitatively new feature, occurring after the initial skin ‘bump’, as comparedto the 2-transition and the TLS results. This feature could be indicative of hemiradialflow, which occurs when a boundary is detected close to the wellbore. Note that theTLS method has smoothed over this feature.In terms of computational cost, there is a significant difference between the TLS andour method, with the former taking usually a few minutes and the latter a few hours,depending on the dataset. However, this trade-off comes with some important inferentialadvantages; namely the uncertainty quantification, the exclusion of non-physical results,and the access to meaningful parameters.In Table 2, we tabulate the AIC (Akaike, 1998), BIC (Schwarz et al., 1978), DIC(Spiegelhalter et al., 2002), and marginal likelihood (Chib and Jeliazkov, 2001) valuesfor all four models. All criteria agree in preferring the 3-transition model, although the4-transition models is a close second. The results for the 3- and 4-transition resultscorroborate the idea that the 4-transition model is over-parameterized, and that threetransitions are sufficient to explain the data in this case.It is worth noting that alternative algorithms, such as the reversible jump MCMC inGreen (1995) can potentially facilitate model selection, by overcoming the requirementof running multiple chains and the subsequent step of using information criteria andBayes factors.
6. Concluding remarks
We have presented a new method for the deconvolution of reservoir well test data, usinga parameterized, physical model for the response function, embedded in a fully Bayesiancontext in order to deal with observational error and other uncertainties. We applied the T. Botsas, I. H. Jermyn and J. A. Cumming
Criterion 1 2 3 4AIC 12781.58 10559.51 T and P .Next steps include developing faster methods to compute the response and, moresubstantially, the extension to the multi-well case (Cumming et al., 2013b).
7. Acknowledgements
This work was performed as part of several Imperial College Joint Industrial Projects(JIP) sponsored by Anadarko, BHP Billiton, BP, Chevron, ConocoPhillips DEA Group, reprint Engie, ENI, Inpex, Nexen CNOOC, Occidental Petroleum, Paradigm, Petro SA, Petro-bras, Saudi Aramco, Schlumberger, Shell, Weatherford and Total.
References
Acosta, L., Ambastha, A. et al. (1994) Thermal well test analysis using an analyticalmulti-region composite reservoir model. In
SPE Annual Technical Conference andExhibition . Society of Petroleum Engineers.Akaike, H. (1998) Information theory and an extension of the maximum likelihood prin-ciple. In
Selected papers of hirotugu akaike , 199–213. Springer.Bourdet, D. (2002)
Well test analysis: the use of advanced interpretation models , vol. 3.Elsevier.Bourdet, D., Whittle, T., Douglas, A. and Pirard, Y. (1983) A new set of type curvessimplifies well test analysis.
World oil , , 95–106.ter Braak, C. J. (2006) A Markov chain Monte Carlo version of the genetic algorithmdifferential evolution: easy Bayesian computing for real parameter spaces. Statisticsand Computing , , 239–249.ter Braak, C. J. and Vrugt, J. A. (2008) Differential evolution Markov chain with snookerupdater and fewer chains. Statistics and Computing , , 435–446.Chib, S. and Jeliazkov, I. (2001) Marginal likelihood from the Metropolis-Hastings out-put. J. Am. Statist. Ass. , , 270–281.Cumming, J., Wooff, D., Whittle, T., Crossman, R., Gringarten, A. et al. (2013a)Assessing the non-uniqueness of the well test interpretation model using deconvolution.In EAGE Annual Conference & Exhibition incorporating SPE Europec . Society ofPetroleum Engineers.Cumming, J., Wooff, D., Whittle, T., Gringarten, A. et al. (2013b) Multiple well decon-volution. In
SPE Annual Technical Conference and Exhibition . Society of PetroleumEngineers. T. Botsas, I. H. Jermyn and J. A. Cumming
Green, P. J. (1995) Reversible jump Markov chain Monte Carlo computation andBayesian model determination.
Biometrika , , 711–732.Gringarten, A. C., Bourdet, D. P., Landel, P. A., Kniazeff, V. J. et al. (1979) A compar-ison between different skin and wellbore storage type-curves for early-time transientanalysis. In SPE Annual Technical Conference and Exhibition . Society of PetroleumEngineers.Gringarten, A. C. et al. (2008) From straight lines to deconvolution: The evolution ofthe state of the art in well test analysis.
SPE Reservoir Evaluation & Engineering , , 41–62.Hartig, F., Minunno, F. and Paul, S. (2019) BayesianTools: General-PurposeMCMC and SMC Samplers and Tools for Bayesian Statistics . URL: https://CRAN.R-project.org/package=BayesianTools . R package version 0.1.6.Horner, D. et al. (1951) Pressure build-up in wells. In .World Petroleum Congress.Satman, A., Eggenschwiler, M., Ramey Jr, H. J. et al. (1980) Interpretation of injec-tion well pressure transient data in thermal oil recovery. In
SPE California RegionalMeeting . Society of Petroleum Engineers.von Schroeter, T., Hollaender, F., Gringarten, A. et al. (2004) Deconvolution of well-testdata as a nonlinear total least-squares problem.
SPE Journal , , 375–390.Schwarz, G. et al. (1978) Estimating the dimension of a model. Ann. Statist. , , 461–464.Spiegelhalter, D. J., Best, N. G., Carlin, B. P. and Van Der Linde, A. (2002) Bayesianmeasures of model complexity and fit. J. R. Statist. Soc. B , , 583–639.Storn, R. and Price, K. (1997) Differential evolution–a simple and efficient heuristicfor global optimization over continuous spaces. Journal of global optimization , ,341–359.Zhang, L., Guo, J. and Liu, Q. (2010) A well test model for stress-sensitive and hetero-geneous reservoirs with non-uniform thicknesses. Petroleum Science ,7