Multi-scale Modeling of 1-D Permeability Fields - Duke University
Multi-scale Modeling of 1-D Permeability Fields Marco A. R. Ferreira1,2 , Zhuoxin Bi3 , Mike West1 , Herbert K. H. Lee1 and David Higdon1,4 1
ISDS - Duke University,
2
DME - UFRJ,
3
Petrotel Inc., 4 Los Alamos National Laboratory.
1. Introduction Permeability plays an important role in subsurface fluid flow studies, being one of the most important quantities for the prediction of fluid flow patterns. The estimation of permeability fields is therefore critical and necessary for the prediction of the behavior of contaminant plumes in aquifers and the production of petroleum from oil fields. In the particular case of production of petroleum from mature fields, part of the available information for the estimation of permeability fields is the production data. To incorporate such information in formal statistical analysis, corresponding likelihood functions for the high-dimensional random field parameters representing the permeability field can be computed with the help of a fluid flow simulator (FFS). Additional information about the permeability field is usually available at different scales of resolution as a result of studies of the geological characteristics of the oil field, well tests, and laboratory measurements. In this paper, in order to incorporate the information available at the different scales of resolution, we use the multi-scale time series model introduced in Ferreira et al. (2001) as a prior for 1-D permeability fields. Estimation of the multi-scale permeability field is then performed using an MCMC algorithm with an embedded FFS running at different scales to incorporate the information given by the production data. We study with simulated data the performance of the proposed approach with respect to the recovery of the original permeability field. 2. General Background A basic concept in subsurface hydrology and petroleum engineering is Darcy’s law which states that in one-dimension, K ∂p v=− , µ ∂x where K = K(x) is the permeability that is used to characterize the porous media, µ is the viscosity of the fluid, p = p(x, t) is the fluid pressure that can be observed at some locations or wells and v = v(x, t) is called Darcy’s velocity or superficial velocity of the fluid. Multiplying Darcy’s velocity by the cross-sectional area A, we obtain the flow rate q, q(x, t) = −
AK ∂p . µ ∂x
Assuming slightly compressible flow, the law of conservation of fluid mass gives the variation of the pressure through time and space, i.e., ∂p ∂ K ∂p ( ) = ct φ , ∂x µ ∂x ∂t
2
Marco A. R. Ferreira et al.
where ct is the total compressibility of the porous medium and φ = φ(x) is the medium porosity which, like permeability, is another parameter characterizing the medium. In order to completely specify the physical system describing the fluid flow, we need initial and boundary conditions. Here, we specify the following three conditions: • p(x, t = 0) = pi , i.e., the initial pressure of the whole field is constant with value pi at time t = 0; •
∂p(x,t) ∂x |x=0
•
AK ∂p(x,t) µ ∂x |x=xwell
= 0, i.e., there is no production at x = 0; = Q, the flow rate at x = xwell is constant equal to Q.
In general, this fluid flow system can not be solved analytically. Solutions are based on numerical approximations such as finite difference methods, generally by partitioning the field into discrete gridblocks. Computer programs that use these numerical approximations to solve this system are called fluid flow simulators (FFS). Depending on the production regime, the output of FFS’s can be either flow rate or pressure at the production well for a given time period that can be up to decades in the future. For references on fluid flow simulators, see King and Datta-Gupta (1998) and Vasco and Datta-Gupta (1999). In this simulation process, the most important parameter used to describe the porous medium is the permeability field which varies by several orders of magnitude more than the porosity field so that φ is tipically treated as a constant. In order to incorporate the production data in the estimation of the permeability field, we must solve the inverse problem, that is, we must estimate the permeability field by historymatching the production data actually observed from the medium. The time- consuming FFS is a must in this process. A review of the inverse problem can be found in Yeh (1986). An early reference on Bayesian methods for the hydrology inverse problem is Neuman and Yakowitz (1979). Floris et al. (1999) discuss the importance of quantifying uncertainty of production forecasts. Craig et al. (1996), Hegstad and Omre (1997), Oliver (1994), Oliver et al. (1997), Vasco and Datta-Gupta (1999) also proposed approaches for the inverse problem that took into account uncertainty aspects. Usually, the inverse problem is ill-posed, i.e., there is never enough data to estimate the permeability field uniquely. In order to find a solution for the inverse problem, it is necessary to impose some regularity conditions on the permeability field. Many investigators in the petroleum industry assume that the permeability field is a Gaussian field process (Cunha et al. 1996, Oliver et al. 1997) or a Markov Random Field (MRF) (Lee et al. 2002). If the permeability field is assumed to be an MRF, then it is possible to solve the inverse problem by using an MCMC algorithm to simulate plausible permeability fields (Lee et al. 2002). At each Metropolis-Hastings step, a permeability field proposal is generated, and the FFS is used to evaluate the likelihood of the proposal in order to decide whether it is accepted or not. See Gamerman (1997) and Gilks et al. (1996) for references on MCMC methodology. The fluid flow simulator only provides an approximate solution. Its accuracy depends on the spatial discretization or the number of gridblocks. As the number of gridblocks grows, the simulator becomes more accurate but slower; as a result, there is a practical upperbound for the number of gridblocks when we are using the fluid flow simulator embedded in an MCMC algorithm. Another important aspect is that the static measurement data on the permeability field may be available at different scales of resolution. The methods proposed so far do not
Multi-Scale permeability modeling
3
address the incorporation of all available information at different scales of resolution in a probabilistically coherent way. Here we propose the use of the Multi-scale Markov Time series models introduced in Ferreira et al. (2001) as a way of coherently incorporating all available information at different scales. It is important to note that the class of Multiscale Markov Time series models induces a different type of regularization than the class of Markov Random Field models. More specifically, this class of multi-scale models induces a long-range type of behavior that can be very useful in the modeling of permeability fields. 3. Multi-scale modeling In general reservoir modeling problems, there are two types of data on the permeabilities: measurement or static data and production data. The measurement information about the permeability field has in general three different levels of aggregation or resolution. The measurement information comes from geological characteristics, well tests and laboratory measurements. The geological data is obtained through seismic tests and gives information on the coarsest scale. The well test data is obtained by producing in the well for a very short period of time and gives information on the intermediate scale. Finally, laboratory measurements are obtained from core samples of the porous media and provide information on the finest scale. We can think about these three sources of information as measurements of the permeability field at the corresponding scales of resolution. The production data are the historic data of oil or water production from the reservoir. Compared with the measurement data, we treat the production data differently because this type of data gives information about the permeability field at several different levels of resolution at the same time. Moreover, the relationship between the production data and the permeability field is highly non-linear. Figure 1 represents graphically the mathematical formulation of our approach to estimating the permeability field. Conceptually, we assume that there is a true continuous permeability field denoted by Xtrue . We observe the production data or curve pobs that is related to the true permeability field through nature laws. In addition, we assume that we can define coarser discretized versions of Xtrue denoted by X0 = coarsest, X1 = intermediate and X2 = finest. Other pieces of information are the measurement data d0 , d1 and d2 on these discretized versions of the true permeability field. Finally, we assume that we have a physical model that approximates well the nature laws and it is implemented in a computer code (FFS) taking as input a discretized permeability field and providing as output the expected production curve. Therefore, the main interest is the incorporation of all available pieces of information in the inference process for the permeability field. We use here the multi-scale time series model introduced in Ferreira et al. (2001) as a prior for the multi-scale permeability field: p(X2 |X1 , µ, φ2 , σ22 , τ2 )p(X1 |X0 , µ, φ1 , σ12 , τ1 )p(X0 |µ, φ0 , σ02 ), where
p(X0 |µ, φ0 , σ02 ) = N (X0 |µ~1n0 , Q0 ) N (X1 |µ~1n1 , V1 )N (X0 |A1 X1 , τ1 In0 ) p(X1 |X0 , µ, φ1 , σ12 , τ1 ) = , N (X0 |µ~1n0 , A1 V1 A01 + τ1 In0 )
and p(X2 |X1 , µ, φ2 , σ22 , τ2 ) =
N (X2 |µ~1n2 , V2 )N (X1 |A2 X2 , τ2 In1 ) , N (X1 |µ~1n1 , A2 V2 A02 + τ2 In1 )
(1)
4
Marco A. R. Ferreira et al.
X
True
X
0
X
1
X
2
Nature FFS p (t) Obs
p (t) 2
Fig. 1. Multi-scale permeability model
d
0
d
1
d
2
Multi-Scale permeability modeling |i−j| where (Q0 )ij = σ02 φ0 /(1 − φ20 ), ~ field and A1 = m−1 1 In1 ⊗ 1m1 is the
5
|i−j| σ12 φ1 /(1
(V1 )ij = − φ21 ), µ is the mean level of the matrix that coarsens the intermediate level X1 by non|i−j| overlapping moving averages. Analogously, (V2 )ij = σ22 φ2 /(1 − φ22 ) and A2 = m−1 2 I n2 ⊗ ~1m2 . We denote the available information by: pobs = observed production data; d0 = measurement of permeability at the coarsest scale; d1 = measurement of permeability at the intermediate scale; d2 = measurement of permeability at the finest scale. We assume that pobs |X1 ∼ N (f (X1 ), σ21 I), d0 |X0 ∼ N (X0 , S0 ), d1 |X1 ∼ N (X1 , S1 ), d2 |X2 ∼ N (X2 , S2 ), where S0 , S1 and S2 are the known covariance matrices of the measurement errors at the respective levels of resolution and f (x) is the production curve computed by the FFS for the permeability field x. Tipically, S0 , S1 and S2 are diagonal matrices or matrices resulting from some geostatistical model applied to the static data. Therefore, the joint distribution of all random quantities is: p(pobs , d0 , d1 , d2 , X0 , X1 , X2 , µ, φ0 , σ02 , φ1 , σ12 , τ1 , φ2 , σ22 , τ2 , σ2 ) ∝ p(pobs |X2 , σ2 )p(d2 |X2 )p(d1 |X1 )p(d0 |X0 )p(X2 |X1 , µ, φ2 , σ22 , τ2 ) p(X1 |X0 , µ, φ1 , σ12 , τ1 ) p(X0 |µ, φ0 , σ02 )p(µ)p(φ2 , σ22 , τ2 )p(φ1 , σ12 , τ1 )p(φ0 , σ02 )p(σ2 ). We want to obtain a summary of the posterior distribution of the permeability field at the different scales: p(X0 , X1 , X2 |pobs , d0 , d1 , d2 ).
(2)
Ultimately, we would really want p(Xtrue |pobs , d0 , d1 , d2 ) ≈ p(X2 |pobs , d0 , d1 , d2 ). Here, we use MCMC to generate a sample from the posterior distribution (2). Each iteration of the MCMC is performed in a cascade way. First, we send all the information about the measurements d0 , d1 and d2 from the finest to the coarsest level. Then, we go back incorporating the production information from the coarsest to the finest level. The method proposed has the following characteristics: • The incorporation of the measurements d0 and d1 from the finest to the coarsest level can be done analytically; • The simulation at coarser levels converges much faster than at finer levels; • The results at coarser levels guide the simulation at finer levels, and the main consequence is faster convergence at finer levels as compared to convergence of traditional MRF schemes. 4. Propagating measurement information from finer to coarser levels In this section we describe how to send the measurement information from finer to coarser levels. Here, we consider only two levels because the generalization to more levels is straightforward. In short, our objective in this section is to obtain p(X0 |d0 , d1 ). For notation simplicity, we omit the dependence of the density functions on the hyperparameters φ0 , σ02 , µ, φ1 , σ12 , τ1 and σ21 . Hence: Z p(X0 |d0 , d1 ) = p(X0 , X1 |d0 , d1 )dX1
6
Marco A. R. Ferreira et al.
∝ = =
∝
Z
Z
p(d0 , d1 |X0 , X1 )p(X0 , X1 )dX1 p(d0 |X0 )p(d1 |X1 )p(X0 )p(X1 |X0 )dX1
N (d0 |X0 , S0 )N (X0 |µ~1n0 , Q0 ) N (X0 |µ~1n0 , A1 V1 A01 + τ1 In0 ) Z × N (d1 |X1 , S1 )N (X0 |A1 X1 , τ1 In0 )N (X1 |µ~1n1 , V1 )dX1 N (d0 |X0 , S0 )N (X0 |µ, Q0 ) N (X0 |b0 , B0 ) N (X0 |A1 µ, A1 V1 A01 + τ1 I)
where B0 = A1 (S1−1 + V1−1 )−1 A01 + τ1 I and b0 = A1 (S1−1 + V1−1 )−1 (S1−1 d1 + V1−1 µ). 5. Incorporating production data from coarsest to finest level In this section we describe how to incorporate production data from the coarsest to the finest level, after having incorporated the measurement data. The main idea of our MCMC algorithm is the following: we generate the coarsest level X0 from p(X0 |pobs , d0 , d1 , d2 ), that is, we integrate out X1 and X2 before generating X0 . Then, we generate X1 from p(X1 |X0 , pobs , d0 , d1 , d2 ) with X2 integrated out. Finally, we generate X2 from its full conditional p(X2 |X1 , X0 , pobs , d0 , d1 , d2 ). In this section, we obtain the expressions of these distributions. The first step of our MCMC algorithm is the simulation of X0 , conditional on the pressure and the measurement data. We assume that the production pobs is independent of the measurements at the finest and intermediate levels given the coarsest level. This is obviously not true, but it is usually a good approximation. The main objective of this assumption is the use of the FFS in the generation of each level of resolution. Let us now obtain the conditional distribution of X0 given pobs , d0 , d1 and d2 : Z Z p(X0 |pobs , d0 , d1 , d2 ) ∝ p(pobs |X2 )p(d2 |X2 )p(d1 |X1 )p(d0 |X0 )p(X2 |X1 ) p(X1 |X0 )p(X0 )dX1 dX2 = p(pobs , d0 , d1 , d2 |X0 )p(X0 ) ∝ p(pobs |X0 )p(X0 |d0 , d1 , d2 ).
(3)
Note that p(X0 |d0 , d1 , d2 ) can be easily obtained by the method outlined in section 4. In addition, a good approximation for p(pobs |X0 ) can be obtained by running the FFS 2 at the coarsest level, that is, p(pobs |X0 ) ≈ N (p R obs R |f (X0 ), σ0 I). In this way, we avoid the computation of the integral in p(pobs |X0 ) = p(pobs |X2 ) p(X2 |X1 )p(X1 |X0 )dX1 dX2 which is very complicated due to the nonlinear relationship between pobs and X2 . In order to obtain the conditional distribution of X1 given X0 , pobs , d0 , d1 and d2 , we assume that pobs is independent of d2 given X1 . Then, p(X1 |X0 , pobs , d0 , d1 , d2 )
∝ p(pobs |X1 , X0 )p(d0 , d1 , d2 |X1 , X0 )p(X1 |X0 ) ∝ p(pobs |X1 )p(X1 |X0 , d0 , d1 , d2 ).
(4)
Multi-Scale permeability modeling
7
But: p(X1 |X0 , d0 , d1 , d2 ) ∝
Z
p(d2 |X2 )p(d1 |X1 )p(d0 |X0 )p(X2 |X1 )p(X1 |X0 )p(X0 )dX2 Z ∝ p(X1 |X0 )p(d1 |X1 ) p(d2 |X2 )p(X2 |X1 )dX2 . (5)
Therefore, p(X1 |X0 , pobs , d0 , d1 , d2 )
∝ p(pobs |X1 )p(X1 |X0 )p(d1 |X1 )
Z
p(d2 |X2 )p(X2 |X1 )dX2 .
R Note that p(d2 |X2 )p(X2 |X1 )dX2 is the operation described in section 4 that brings the measurement information about the finest level to the intermediate level. In addition, p(pobs |X1 ) can be approximated by running the FFS at the intermediate level and assuming that p(pobs |X1 ) ≈ N (pobs |f (X1 ), σ21 I).
Porous media Production well
1280 m
Fig. 2. Subsurface porous media.
Finally, the conditional distribution of X2 given X1 , X0 , pobs , d0 , d1 and d2 is p(X2 |X1 , X0 , pobs , d0 , d1 , d2 )
∝ p(pobs |X2 )p(X2 |X1 , X0 , d0 , d1 , d2 ) ∝ p(pobs |X2 )p(X2 |X1 )p(d2 |X2 ).
(6)
6. Application In this section we present an application with a synthetic 1-D permeability field. The design of the field is presented in Figure 2. The length of the field is equal to 1280 meters and there is only one production well located in the right end of the field. Here, we assume that the finest discretized level X2 is equal to Xtrue and it is unfeasible to run the FFS at the finest level. Therefore, it is not possible to perform inference on X2 . In this application, we focus on the inference process on the coarser versions X0 and X1 . In subsection 6.1 we discuss the specification of the prior distribution. In subsection 6.2 we present the setup of the application. The results of the analysis are discussed in subsection 6.3.
8
Marco A. R. Ferreira et al.
6.1. Prior specification In this subsection we discuss the specification of the prior distribution for the parameters of the model. In order to reduce the dimension of the prior specification problem, we assume a priori independent knowledge about the parameters, allowing us to focus on the specification of the prior for each parameter individually. There are several aspects peculiar to our approach to the 1-D fluid flow problem that allow us to set informative priors for several parameters of interest. The most important of these characteristics is that the permeability field is modeled on the logarithm scale. Hence the variances σ02 , σ12 and τ1 are related to the amount of relative variability that we are willing to allow. The autoregressive coefficients φ0 and φ1 are related to the amount of smoothness that we want to impose at each level of the log-permeability field. The overall log-permeability mean µ is assigned a prior with mass in the region where the logpermeability makes physical sense. The variances σ0 and σ1 are assigned vague priors. More specifically, the priors are φ0 ∼ N (0.0, 0.5), σ02 ∼ IG(3.0, 1.1), φ1 ∼ N (0.9, 0.01), σ12 ∼ IG(57, 0.14) µ ∼ N (5.0, 1.21), σ0 ∼ IG(5 × 10−5 , 5 × 10−5 ), σ1 ∼ IG(5 × 10−5 , 5 × 10−5 ), τ1 ∼ IG(57, 0.14). 6.2. Problem Setup
(a)
(b)
(c) Fig. 3. Original field and coarser versions. (a) Original field - 512 gridblocks; (b) Intermediate field 64 gridblocks; (c) Coarsest field - 8 gridblocks.
Figure 3 presents the original synthetic permeability field with 512 gridblocks and coarser versions with 64 and 8 gridblocks. The i-th gridblock of the synthetic permeability field was set equal to exp(−3.1710−5i2 + 0.01625i + 3.91), i = 0, . . . , 511. The other characteristics of the porous media are assumed constant and known. More specifically, in this example the following values were used: • Single phase;
ν = Viscosity = 1.0cp;
• pin = Initial pressure = 1500psi; • φ = Porosity = 0.25;
Multi-Scale permeability modeling Table 1. Measurement data at intermediate level of resolution j 1 9 17 25 33 41 49 57 d1j 3.97 4.86 5.50 5.88 5.99 5.85 5.45 4.78
9
64 3.99
• Q = Production = 2rb/day; • Cross sectional area = 100m2 ;
1400 1300
1350
pressure
1450
1500
• Compressibility = 0.4 × 10−5 psi−1 .
0.0
0.5
1.0
1.5
2.0
time
Fig. 4. Production data – Pressure curve through time.
The production data in this application are given by the pressure curve presented in Figure 6.2 for the first two days of production, with variance of the measures of pressure equal to σ2 = 1.0. The measurements at the coarse level are d0 = (4.4, 5.1, 5.6, 5.9, 6.0, 5.7, 5.2, 4.3) and S0 = 0.12 I. Table 1 presents independent measurements at the intermediate level observed only at 9 gridblocks with variance equal to 0.004.
6.3. Results In this subsection we discuss the results of the application of the methodology proposed to perform inference about the permeability field through the analysis of the production and measurement data. Table 2 presents posterior summaries for the hyperparameters of the multi-scale model. The estimates of φ0 , σ02 , φ1 and σ12 are evidence of the higher degree of smoothing/regularization obtained at the intermediate level when compared with the regularization at the coarsest level. The estimate for τ1 shows that each coarsest gridblock is very close to the average of the corresponding intermediate gridblocks. Finally, the σ20 and σ21 evidence that history matching is more precise at the intermediate level than at the coarsest level. Moreover, as the original variance of the measures of the precision curver was σ2 = 1.0, there is a fair amount of overfitting of the production curve in both levels of resolution.
10
Marco A. R. Ferreira et al. Table 2. Posterior summaries for the hyperparameters Mean Standard deviation φ0 0.4483 0.3097 σ02 0.4416 0.2207 φ1 0.8630 0.0319 σ12 0.0024 0.0005 τ1 0.0019 0.0003 σ20 0.9536 0.0874 σ21 0.8915 0.0801 µ 5.0680 0.4810
Figure 5 presents the coarsest version of the original field, posterior mean and Gibbs sample for the coarse level of resolution. Figure 5c represents the Gibbs sample, each column corresponding to one gridblock. In addition, each line in Figure 5c represents a realization of the log-permeability field, the bottom line representing the initial log-permeability field set constant equal to 5 log milidarcies, and the top line representing the last realization. As we can see in Figure 5c, the Gibbs sampling converges quite rapidly in less than 1000 iterations. Figure 5a represents the coarsest version of the original field. Figure 5b is the posterior mean of the coarsest level obtained by averaging the Gibbs sample realizations after a burnin period of 5000 iterations. The posterior mean is very close to the coarsest level of the original permeability field indicating that our methodology is working well. Analogously, Figure 6 presents the intermediate version of the original field, posterior mean and Gibbs sample for the intermediate level of resolution. The conclusions are similar to those obtained for the coarsest level, the main difference being that the estimation at the intermediate level of resolution recovers a higher degree of detail than the estimation at the coarsest level. As stated by Oliver et al. (1997), the main advantage of the Bayesian framework coupled with MCMC technology is the characterization of the uncertainty present in the inference about the permeability field. A useful summary of this uncertainty is the posterior variance of each gridblock. Figure 7 presents the posterior variance across the field at the coarsest and intermediate levels of resolution. In Figure 7a, the solid line represents the posterior variance at the coarsest level while the dashed line represents the variance of the measurements d 0i . The difference between the solid and the dashed lines is due to the use of the multi-scale 1-D model as prior for the log-permeability field and to the incorporation of the production data in the inference process. As the production well is located in the right end of the field, the incorporation of the production data results in smaller posterior variances for the gridblocks located on the right of the field. In particular, the ratio between the posterior variances of gridblocks 1 (left end) and 8 (right end) is equal to 50, giving evidence about the reduction in uncertainty due to the proximity of the production well. In figure 7b, the solid line represents the posterior variance at the intermediate level. The dashed line is the variance of the measurement data obtained for gridblocks 1, 9, 17, 25, 33, 41, 49, 57 and 64. Note that the posterior variances at these gridblocks are reasonably smaller than the measurement variance. Moreover, although the posterior variances of the gridblocks without measurement information are in general higher than the posterior variances of the gridblocks with measurement information, the multi-scale model transfers information between gridblocks and levels of resolution. As the measurement information
Multi-Scale permeability modeling
11
(a)
10000 0
5000
Iterations
15000
20000
(b)
2
4
6
8
Gridblocks
(c) Fig. 5. Coarsest level of resolution. (a) Version of original field; (b) Posterior mean; (c) Gibbs Sample.
is equally spaced, the transfer of information along the field causes periodic behavior of the posterior variances. Finally, as in the coarsest level of resolution, the posterior variances of the gridblocks on the right end of the field are smaller because of the use of the production data.
7. Conclusions and discussion In this paper, we used a 1-D multi-scale model as a prior for a permeability field. We developed an MCMC algorithm embedded in the Bayesian framework such that a Fluid Flow Simulator can be used at each level of resolution to incorporate the production information. Moreover, the algorithm is performed in a cascade way, from coarser to finer levels. In this way, the realizations at coarser levels guide the realizations at the finer levels leading to a more efficient exploration of the posterior distribution of the multi-scale permeability field. The final result is a multi-scale framework that incorporates data at different resolutions and/or of different types in a coherent and efficient way. We applied the proposed methodology to synthetic data. As expected from intuition, production data gives more information on permeabilities closer to the production well. Finally, the use of data from several different sources allows a very good recovery of the original permeability field. Due to the restriction to 1 dimension, applications to real world problems do not abound.
12
Marco A. R. Ferreira et al.
(a)
10000 0
5000
Iterations
15000
20000
(b)
0
10
20
30
40
50
60
Gridblocks
(c) Fig. 6. Intermediate level of resolution. (a) Version of original field; (b) Posterior mean; (c) Gibbs Sample.
One example of potential application of the proposed methodology is to petroleum fields that exist in the bed of rivers. But the main ideas of the paper survive: the use of a multiscale model defined from coarser to finer levels as prior for the permeability field and the use of the Bayesian paradigm to perform inference about the parameters of interest. Ongoing research is the extension of this work to the 2 dimensional case. The main questions are the development of theoretical multi-scale Markov Random Fields and the development of estimation and simulation procedures. A second step is the adjustment of these theoretical developments to hydrology applications. A. MCMC details The estimation of the parameters of the model is done using an MCMC algorithm. The full conditionals for φ0 and σ02 are trivially obtained as discussed in Chib and Greenberg (1994) and McCulloch and Tsay (1994). The full conditionals for φ1 , σ12 , τ1 , φ2 , σ22 and τ2 are as discussed in Ferreira et al. (2001). Within the MCMC scheme, the various resolution levels of the permeability field are updated with block Metropolis-Hastings steps. In addition, we integrate out the finer levels of resolution when generating a given level of the permeability field. Therefore, the generation of the current level i depends only on the immediate coarser level of resolution
13
0.0
0.002
0.004
0.006
0.008
0.010
Multi-Scale permeability modeling
1
2
3
4
5
6
7
8
0.0015
0.0025
0.0035
0.0045
(a)
0
10
20
30
40
50
60
(b) Fig. 7. Posterior variance across the field. (a) coarse level. (b) intermediate level. Solid line = posterior variance. Dashed line = variance of the measurement data
i − 1. The full conditional for σ21 is equal to IG(n∗ , n∗ s∗ /2), where n∗ = n + T and n∗ s∗ = n s + (pobs − f (X1 ))0 (pobs − f (X1 )). The full conditional for σ20 is obtained analogously. References Chib, S. and Greenberg, E. (1994). Bayes inference in regression models with arma(p,q) errors. Journal of Econometrics, 64:183–206. Craig, P. S., Goldstein, M., Seheult, A. H., and Smith, J. A. (1996). Bayes linear strategies for history matching of hydrocarbon reservoirs. In Bernardo, J. M., Berger, J. O., Dawid, A. P., and Smith, A. F. M., editors, Bayesian Statistics 5, pages 69–95. Clarendon Press, Oxford. Cunha, L. B., Oliver, D. S., Redner, R. A., and Reynolds, A. C. (1996). A hybrid markov chain monte carlo method for generating permeability fields conditioned to multiwell
14
Marco A. R. Ferreira et al.
pressure data and prior information. In SPE Annual Technical Conference and Exhibition, Denver. SPE. Ferreira, M. A. R., West, M., Lee, H., and Higdon, D. (2001). A class of multi-scale time series models. Technical report, Duke University - ISDS. Floris, F. J. T., Bush, M. D., Cuypers, M., Roggero, F., and Syversveen, A.-R. (1999). Comparison of production forecast uncertainty quantification methods — an integrated study. 1st Symposium on Petroleum Geostatistics, Toulouse, 20–23 April 1999. Gamerman, D. (1997). Markov Chain Monte Carlo. Chapman Hall. Gilks, W. R., Richardson, S., and Spiegelhalter, D. J., editors (1996). Markov Chain Monte Carlo in Practice. Chapman Hall. Hegstad, B. K. and Omre, H. (1997). Uncertainty assessment in history matching and forecasting. In E. Y. Baafi and N. A. Schofield, editor, Geostatistics Wollongong ’96, Vol. 1., pages 585–596. Kluwer Academic Publishers. King, M. J. and Datta-Gupta, A. (1998). Streamline simulation: a current perspective. In Situ, 22(1):91–140. Lee, H., Higdon, D., Bi, Z., Ferreira, M. A. R., and West, M. (2002). Markov random field models for high-dimensional parameters in simulations of fluid flow in porous media. Technometrics. To appear. McCulloch, R. E. and Tsay, R. S. (1994). Bayesian analysis of autoregressive time series via the gibbs sampler. Journal of time series analysis, 15:235–250. Neuman, S. P. and Yakowitz, S. (1979). A statistical approach to the problem of aquifer hydrology: 1. theory. Water Resources Research, 15:845–860. Oliver, D. S. (1994). Incorporation of transient pressure data into reservoir characterization. In Situ, 18:243–275. Oliver, D. S., Cunha, L. B., and Reynolds, A. C. (1997). Markov chain monte carlo methods for conditioning a permeability field to pressure data. Mathematical Geology, 29:61–91. Vasco, D. W. and Datta-Gupta, A. (1999). Asymptotic solutions for solute transport: A formalism for tracer tomography. Water Resources Research, 35(1):1–16. Yeh, W. W. (1986). Review of parameter identification in groundwater hydrology: the inverse problem. Water Resources Research, 22:95–108.
ISDS - Duke University,
2
DME - UFRJ,
3
Petrotel Inc., 4 Los Alamos National Laboratory.
1. Introduction Permeability plays an important role in subsurface fluid flow studies, being one of the most important quantities for the prediction of fluid flow patterns. The estimation of permeability fields is therefore critical and necessary for the prediction of the behavior of contaminant plumes in aquifers and the production of petroleum from oil fields. In the particular case of production of petroleum from mature fields, part of the available information for the estimation of permeability fields is the production data. To incorporate such information in formal statistical analysis, corresponding likelihood functions for the high-dimensional random field parameters representing the permeability field can be computed with the help of a fluid flow simulator (FFS). Additional information about the permeability field is usually available at different scales of resolution as a result of studies of the geological characteristics of the oil field, well tests, and laboratory measurements. In this paper, in order to incorporate the information available at the different scales of resolution, we use the multi-scale time series model introduced in Ferreira et al. (2001) as a prior for 1-D permeability fields. Estimation of the multi-scale permeability field is then performed using an MCMC algorithm with an embedded FFS running at different scales to incorporate the information given by the production data. We study with simulated data the performance of the proposed approach with respect to the recovery of the original permeability field. 2. General Background A basic concept in subsurface hydrology and petroleum engineering is Darcy’s law which states that in one-dimension, K ∂p v=− , µ ∂x where K = K(x) is the permeability that is used to characterize the porous media, µ is the viscosity of the fluid, p = p(x, t) is the fluid pressure that can be observed at some locations or wells and v = v(x, t) is called Darcy’s velocity or superficial velocity of the fluid. Multiplying Darcy’s velocity by the cross-sectional area A, we obtain the flow rate q, q(x, t) = −
AK ∂p . µ ∂x
Assuming slightly compressible flow, the law of conservation of fluid mass gives the variation of the pressure through time and space, i.e., ∂p ∂ K ∂p ( ) = ct φ , ∂x µ ∂x ∂t
2
Marco A. R. Ferreira et al.
where ct is the total compressibility of the porous medium and φ = φ(x) is the medium porosity which, like permeability, is another parameter characterizing the medium. In order to completely specify the physical system describing the fluid flow, we need initial and boundary conditions. Here, we specify the following three conditions: • p(x, t = 0) = pi , i.e., the initial pressure of the whole field is constant with value pi at time t = 0; •
∂p(x,t) ∂x |x=0
•
AK ∂p(x,t) µ ∂x |x=xwell
= 0, i.e., there is no production at x = 0; = Q, the flow rate at x = xwell is constant equal to Q.
In general, this fluid flow system can not be solved analytically. Solutions are based on numerical approximations such as finite difference methods, generally by partitioning the field into discrete gridblocks. Computer programs that use these numerical approximations to solve this system are called fluid flow simulators (FFS). Depending on the production regime, the output of FFS’s can be either flow rate or pressure at the production well for a given time period that can be up to decades in the future. For references on fluid flow simulators, see King and Datta-Gupta (1998) and Vasco and Datta-Gupta (1999). In this simulation process, the most important parameter used to describe the porous medium is the permeability field which varies by several orders of magnitude more than the porosity field so that φ is tipically treated as a constant. In order to incorporate the production data in the estimation of the permeability field, we must solve the inverse problem, that is, we must estimate the permeability field by historymatching the production data actually observed from the medium. The time- consuming FFS is a must in this process. A review of the inverse problem can be found in Yeh (1986). An early reference on Bayesian methods for the hydrology inverse problem is Neuman and Yakowitz (1979). Floris et al. (1999) discuss the importance of quantifying uncertainty of production forecasts. Craig et al. (1996), Hegstad and Omre (1997), Oliver (1994), Oliver et al. (1997), Vasco and Datta-Gupta (1999) also proposed approaches for the inverse problem that took into account uncertainty aspects. Usually, the inverse problem is ill-posed, i.e., there is never enough data to estimate the permeability field uniquely. In order to find a solution for the inverse problem, it is necessary to impose some regularity conditions on the permeability field. Many investigators in the petroleum industry assume that the permeability field is a Gaussian field process (Cunha et al. 1996, Oliver et al. 1997) or a Markov Random Field (MRF) (Lee et al. 2002). If the permeability field is assumed to be an MRF, then it is possible to solve the inverse problem by using an MCMC algorithm to simulate plausible permeability fields (Lee et al. 2002). At each Metropolis-Hastings step, a permeability field proposal is generated, and the FFS is used to evaluate the likelihood of the proposal in order to decide whether it is accepted or not. See Gamerman (1997) and Gilks et al. (1996) for references on MCMC methodology. The fluid flow simulator only provides an approximate solution. Its accuracy depends on the spatial discretization or the number of gridblocks. As the number of gridblocks grows, the simulator becomes more accurate but slower; as a result, there is a practical upperbound for the number of gridblocks when we are using the fluid flow simulator embedded in an MCMC algorithm. Another important aspect is that the static measurement data on the permeability field may be available at different scales of resolution. The methods proposed so far do not
Multi-Scale permeability modeling
3
address the incorporation of all available information at different scales of resolution in a probabilistically coherent way. Here we propose the use of the Multi-scale Markov Time series models introduced in Ferreira et al. (2001) as a way of coherently incorporating all available information at different scales. It is important to note that the class of Multiscale Markov Time series models induces a different type of regularization than the class of Markov Random Field models. More specifically, this class of multi-scale models induces a long-range type of behavior that can be very useful in the modeling of permeability fields. 3. Multi-scale modeling In general reservoir modeling problems, there are two types of data on the permeabilities: measurement or static data and production data. The measurement information about the permeability field has in general three different levels of aggregation or resolution. The measurement information comes from geological characteristics, well tests and laboratory measurements. The geological data is obtained through seismic tests and gives information on the coarsest scale. The well test data is obtained by producing in the well for a very short period of time and gives information on the intermediate scale. Finally, laboratory measurements are obtained from core samples of the porous media and provide information on the finest scale. We can think about these three sources of information as measurements of the permeability field at the corresponding scales of resolution. The production data are the historic data of oil or water production from the reservoir. Compared with the measurement data, we treat the production data differently because this type of data gives information about the permeability field at several different levels of resolution at the same time. Moreover, the relationship between the production data and the permeability field is highly non-linear. Figure 1 represents graphically the mathematical formulation of our approach to estimating the permeability field. Conceptually, we assume that there is a true continuous permeability field denoted by Xtrue . We observe the production data or curve pobs that is related to the true permeability field through nature laws. In addition, we assume that we can define coarser discretized versions of Xtrue denoted by X0 = coarsest, X1 = intermediate and X2 = finest. Other pieces of information are the measurement data d0 , d1 and d2 on these discretized versions of the true permeability field. Finally, we assume that we have a physical model that approximates well the nature laws and it is implemented in a computer code (FFS) taking as input a discretized permeability field and providing as output the expected production curve. Therefore, the main interest is the incorporation of all available pieces of information in the inference process for the permeability field. We use here the multi-scale time series model introduced in Ferreira et al. (2001) as a prior for the multi-scale permeability field: p(X2 |X1 , µ, φ2 , σ22 , τ2 )p(X1 |X0 , µ, φ1 , σ12 , τ1 )p(X0 |µ, φ0 , σ02 ), where
p(X0 |µ, φ0 , σ02 ) = N (X0 |µ~1n0 , Q0 ) N (X1 |µ~1n1 , V1 )N (X0 |A1 X1 , τ1 In0 ) p(X1 |X0 , µ, φ1 , σ12 , τ1 ) = , N (X0 |µ~1n0 , A1 V1 A01 + τ1 In0 )
and p(X2 |X1 , µ, φ2 , σ22 , τ2 ) =
N (X2 |µ~1n2 , V2 )N (X1 |A2 X2 , τ2 In1 ) , N (X1 |µ~1n1 , A2 V2 A02 + τ2 In1 )
(1)
4
Marco A. R. Ferreira et al.
X
True
X
0
X
1
X
2
Nature FFS p (t) Obs
p (t) 2
Fig. 1. Multi-scale permeability model
d
0
d
1
d
2
Multi-Scale permeability modeling |i−j| where (Q0 )ij = σ02 φ0 /(1 − φ20 ), ~ field and A1 = m−1 1 In1 ⊗ 1m1 is the
5
|i−j| σ12 φ1 /(1
(V1 )ij = − φ21 ), µ is the mean level of the matrix that coarsens the intermediate level X1 by non|i−j| overlapping moving averages. Analogously, (V2 )ij = σ22 φ2 /(1 − φ22 ) and A2 = m−1 2 I n2 ⊗ ~1m2 . We denote the available information by: pobs = observed production data; d0 = measurement of permeability at the coarsest scale; d1 = measurement of permeability at the intermediate scale; d2 = measurement of permeability at the finest scale. We assume that pobs |X1 ∼ N (f (X1 ), σ21 I), d0 |X0 ∼ N (X0 , S0 ), d1 |X1 ∼ N (X1 , S1 ), d2 |X2 ∼ N (X2 , S2 ), where S0 , S1 and S2 are the known covariance matrices of the measurement errors at the respective levels of resolution and f (x) is the production curve computed by the FFS for the permeability field x. Tipically, S0 , S1 and S2 are diagonal matrices or matrices resulting from some geostatistical model applied to the static data. Therefore, the joint distribution of all random quantities is: p(pobs , d0 , d1 , d2 , X0 , X1 , X2 , µ, φ0 , σ02 , φ1 , σ12 , τ1 , φ2 , σ22 , τ2 , σ2 ) ∝ p(pobs |X2 , σ2 )p(d2 |X2 )p(d1 |X1 )p(d0 |X0 )p(X2 |X1 , µ, φ2 , σ22 , τ2 ) p(X1 |X0 , µ, φ1 , σ12 , τ1 ) p(X0 |µ, φ0 , σ02 )p(µ)p(φ2 , σ22 , τ2 )p(φ1 , σ12 , τ1 )p(φ0 , σ02 )p(σ2 ). We want to obtain a summary of the posterior distribution of the permeability field at the different scales: p(X0 , X1 , X2 |pobs , d0 , d1 , d2 ).
(2)
Ultimately, we would really want p(Xtrue |pobs , d0 , d1 , d2 ) ≈ p(X2 |pobs , d0 , d1 , d2 ). Here, we use MCMC to generate a sample from the posterior distribution (2). Each iteration of the MCMC is performed in a cascade way. First, we send all the information about the measurements d0 , d1 and d2 from the finest to the coarsest level. Then, we go back incorporating the production information from the coarsest to the finest level. The method proposed has the following characteristics: • The incorporation of the measurements d0 and d1 from the finest to the coarsest level can be done analytically; • The simulation at coarser levels converges much faster than at finer levels; • The results at coarser levels guide the simulation at finer levels, and the main consequence is faster convergence at finer levels as compared to convergence of traditional MRF schemes. 4. Propagating measurement information from finer to coarser levels In this section we describe how to send the measurement information from finer to coarser levels. Here, we consider only two levels because the generalization to more levels is straightforward. In short, our objective in this section is to obtain p(X0 |d0 , d1 ). For notation simplicity, we omit the dependence of the density functions on the hyperparameters φ0 , σ02 , µ, φ1 , σ12 , τ1 and σ21 . Hence: Z p(X0 |d0 , d1 ) = p(X0 , X1 |d0 , d1 )dX1
6
Marco A. R. Ferreira et al.
∝ = =
∝
Z
Z
p(d0 , d1 |X0 , X1 )p(X0 , X1 )dX1 p(d0 |X0 )p(d1 |X1 )p(X0 )p(X1 |X0 )dX1
N (d0 |X0 , S0 )N (X0 |µ~1n0 , Q0 ) N (X0 |µ~1n0 , A1 V1 A01 + τ1 In0 ) Z × N (d1 |X1 , S1 )N (X0 |A1 X1 , τ1 In0 )N (X1 |µ~1n1 , V1 )dX1 N (d0 |X0 , S0 )N (X0 |µ, Q0 ) N (X0 |b0 , B0 ) N (X0 |A1 µ, A1 V1 A01 + τ1 I)
where B0 = A1 (S1−1 + V1−1 )−1 A01 + τ1 I and b0 = A1 (S1−1 + V1−1 )−1 (S1−1 d1 + V1−1 µ). 5. Incorporating production data from coarsest to finest level In this section we describe how to incorporate production data from the coarsest to the finest level, after having incorporated the measurement data. The main idea of our MCMC algorithm is the following: we generate the coarsest level X0 from p(X0 |pobs , d0 , d1 , d2 ), that is, we integrate out X1 and X2 before generating X0 . Then, we generate X1 from p(X1 |X0 , pobs , d0 , d1 , d2 ) with X2 integrated out. Finally, we generate X2 from its full conditional p(X2 |X1 , X0 , pobs , d0 , d1 , d2 ). In this section, we obtain the expressions of these distributions. The first step of our MCMC algorithm is the simulation of X0 , conditional on the pressure and the measurement data. We assume that the production pobs is independent of the measurements at the finest and intermediate levels given the coarsest level. This is obviously not true, but it is usually a good approximation. The main objective of this assumption is the use of the FFS in the generation of each level of resolution. Let us now obtain the conditional distribution of X0 given pobs , d0 , d1 and d2 : Z Z p(X0 |pobs , d0 , d1 , d2 ) ∝ p(pobs |X2 )p(d2 |X2 )p(d1 |X1 )p(d0 |X0 )p(X2 |X1 ) p(X1 |X0 )p(X0 )dX1 dX2 = p(pobs , d0 , d1 , d2 |X0 )p(X0 ) ∝ p(pobs |X0 )p(X0 |d0 , d1 , d2 ).
(3)
Note that p(X0 |d0 , d1 , d2 ) can be easily obtained by the method outlined in section 4. In addition, a good approximation for p(pobs |X0 ) can be obtained by running the FFS 2 at the coarsest level, that is, p(pobs |X0 ) ≈ N (p R obs R |f (X0 ), σ0 I). In this way, we avoid the computation of the integral in p(pobs |X0 ) = p(pobs |X2 ) p(X2 |X1 )p(X1 |X0 )dX1 dX2 which is very complicated due to the nonlinear relationship between pobs and X2 . In order to obtain the conditional distribution of X1 given X0 , pobs , d0 , d1 and d2 , we assume that pobs is independent of d2 given X1 . Then, p(X1 |X0 , pobs , d0 , d1 , d2 )
∝ p(pobs |X1 , X0 )p(d0 , d1 , d2 |X1 , X0 )p(X1 |X0 ) ∝ p(pobs |X1 )p(X1 |X0 , d0 , d1 , d2 ).
(4)
Multi-Scale permeability modeling
7
But: p(X1 |X0 , d0 , d1 , d2 ) ∝
Z
p(d2 |X2 )p(d1 |X1 )p(d0 |X0 )p(X2 |X1 )p(X1 |X0 )p(X0 )dX2 Z ∝ p(X1 |X0 )p(d1 |X1 ) p(d2 |X2 )p(X2 |X1 )dX2 . (5)
Therefore, p(X1 |X0 , pobs , d0 , d1 , d2 )
∝ p(pobs |X1 )p(X1 |X0 )p(d1 |X1 )
Z
p(d2 |X2 )p(X2 |X1 )dX2 .
R Note that p(d2 |X2 )p(X2 |X1 )dX2 is the operation described in section 4 that brings the measurement information about the finest level to the intermediate level. In addition, p(pobs |X1 ) can be approximated by running the FFS at the intermediate level and assuming that p(pobs |X1 ) ≈ N (pobs |f (X1 ), σ21 I).
Porous media Production well
1280 m
Fig. 2. Subsurface porous media.
Finally, the conditional distribution of X2 given X1 , X0 , pobs , d0 , d1 and d2 is p(X2 |X1 , X0 , pobs , d0 , d1 , d2 )
∝ p(pobs |X2 )p(X2 |X1 , X0 , d0 , d1 , d2 ) ∝ p(pobs |X2 )p(X2 |X1 )p(d2 |X2 ).
(6)
6. Application In this section we present an application with a synthetic 1-D permeability field. The design of the field is presented in Figure 2. The length of the field is equal to 1280 meters and there is only one production well located in the right end of the field. Here, we assume that the finest discretized level X2 is equal to Xtrue and it is unfeasible to run the FFS at the finest level. Therefore, it is not possible to perform inference on X2 . In this application, we focus on the inference process on the coarser versions X0 and X1 . In subsection 6.1 we discuss the specification of the prior distribution. In subsection 6.2 we present the setup of the application. The results of the analysis are discussed in subsection 6.3.
8
Marco A. R. Ferreira et al.
6.1. Prior specification In this subsection we discuss the specification of the prior distribution for the parameters of the model. In order to reduce the dimension of the prior specification problem, we assume a priori independent knowledge about the parameters, allowing us to focus on the specification of the prior for each parameter individually. There are several aspects peculiar to our approach to the 1-D fluid flow problem that allow us to set informative priors for several parameters of interest. The most important of these characteristics is that the permeability field is modeled on the logarithm scale. Hence the variances σ02 , σ12 and τ1 are related to the amount of relative variability that we are willing to allow. The autoregressive coefficients φ0 and φ1 are related to the amount of smoothness that we want to impose at each level of the log-permeability field. The overall log-permeability mean µ is assigned a prior with mass in the region where the logpermeability makes physical sense. The variances σ0 and σ1 are assigned vague priors. More specifically, the priors are φ0 ∼ N (0.0, 0.5), σ02 ∼ IG(3.0, 1.1), φ1 ∼ N (0.9, 0.01), σ12 ∼ IG(57, 0.14) µ ∼ N (5.0, 1.21), σ0 ∼ IG(5 × 10−5 , 5 × 10−5 ), σ1 ∼ IG(5 × 10−5 , 5 × 10−5 ), τ1 ∼ IG(57, 0.14). 6.2. Problem Setup
(a)
(b)
(c) Fig. 3. Original field and coarser versions. (a) Original field - 512 gridblocks; (b) Intermediate field 64 gridblocks; (c) Coarsest field - 8 gridblocks.
Figure 3 presents the original synthetic permeability field with 512 gridblocks and coarser versions with 64 and 8 gridblocks. The i-th gridblock of the synthetic permeability field was set equal to exp(−3.1710−5i2 + 0.01625i + 3.91), i = 0, . . . , 511. The other characteristics of the porous media are assumed constant and known. More specifically, in this example the following values were used: • Single phase;
ν = Viscosity = 1.0cp;
• pin = Initial pressure = 1500psi; • φ = Porosity = 0.25;
Multi-Scale permeability modeling Table 1. Measurement data at intermediate level of resolution j 1 9 17 25 33 41 49 57 d1j 3.97 4.86 5.50 5.88 5.99 5.85 5.45 4.78
9
64 3.99
• Q = Production = 2rb/day; • Cross sectional area = 100m2 ;
1400 1300
1350
pressure
1450
1500
• Compressibility = 0.4 × 10−5 psi−1 .
0.0
0.5
1.0
1.5
2.0
time
Fig. 4. Production data – Pressure curve through time.
The production data in this application are given by the pressure curve presented in Figure 6.2 for the first two days of production, with variance of the measures of pressure equal to σ2 = 1.0. The measurements at the coarse level are d0 = (4.4, 5.1, 5.6, 5.9, 6.0, 5.7, 5.2, 4.3) and S0 = 0.12 I. Table 1 presents independent measurements at the intermediate level observed only at 9 gridblocks with variance equal to 0.004.
6.3. Results In this subsection we discuss the results of the application of the methodology proposed to perform inference about the permeability field through the analysis of the production and measurement data. Table 2 presents posterior summaries for the hyperparameters of the multi-scale model. The estimates of φ0 , σ02 , φ1 and σ12 are evidence of the higher degree of smoothing/regularization obtained at the intermediate level when compared with the regularization at the coarsest level. The estimate for τ1 shows that each coarsest gridblock is very close to the average of the corresponding intermediate gridblocks. Finally, the σ20 and σ21 evidence that history matching is more precise at the intermediate level than at the coarsest level. Moreover, as the original variance of the measures of the precision curver was σ2 = 1.0, there is a fair amount of overfitting of the production curve in both levels of resolution.
10
Marco A. R. Ferreira et al. Table 2. Posterior summaries for the hyperparameters Mean Standard deviation φ0 0.4483 0.3097 σ02 0.4416 0.2207 φ1 0.8630 0.0319 σ12 0.0024 0.0005 τ1 0.0019 0.0003 σ20 0.9536 0.0874 σ21 0.8915 0.0801 µ 5.0680 0.4810
Figure 5 presents the coarsest version of the original field, posterior mean and Gibbs sample for the coarse level of resolution. Figure 5c represents the Gibbs sample, each column corresponding to one gridblock. In addition, each line in Figure 5c represents a realization of the log-permeability field, the bottom line representing the initial log-permeability field set constant equal to 5 log milidarcies, and the top line representing the last realization. As we can see in Figure 5c, the Gibbs sampling converges quite rapidly in less than 1000 iterations. Figure 5a represents the coarsest version of the original field. Figure 5b is the posterior mean of the coarsest level obtained by averaging the Gibbs sample realizations after a burnin period of 5000 iterations. The posterior mean is very close to the coarsest level of the original permeability field indicating that our methodology is working well. Analogously, Figure 6 presents the intermediate version of the original field, posterior mean and Gibbs sample for the intermediate level of resolution. The conclusions are similar to those obtained for the coarsest level, the main difference being that the estimation at the intermediate level of resolution recovers a higher degree of detail than the estimation at the coarsest level. As stated by Oliver et al. (1997), the main advantage of the Bayesian framework coupled with MCMC technology is the characterization of the uncertainty present in the inference about the permeability field. A useful summary of this uncertainty is the posterior variance of each gridblock. Figure 7 presents the posterior variance across the field at the coarsest and intermediate levels of resolution. In Figure 7a, the solid line represents the posterior variance at the coarsest level while the dashed line represents the variance of the measurements d 0i . The difference between the solid and the dashed lines is due to the use of the multi-scale 1-D model as prior for the log-permeability field and to the incorporation of the production data in the inference process. As the production well is located in the right end of the field, the incorporation of the production data results in smaller posterior variances for the gridblocks located on the right of the field. In particular, the ratio between the posterior variances of gridblocks 1 (left end) and 8 (right end) is equal to 50, giving evidence about the reduction in uncertainty due to the proximity of the production well. In figure 7b, the solid line represents the posterior variance at the intermediate level. The dashed line is the variance of the measurement data obtained for gridblocks 1, 9, 17, 25, 33, 41, 49, 57 and 64. Note that the posterior variances at these gridblocks are reasonably smaller than the measurement variance. Moreover, although the posterior variances of the gridblocks without measurement information are in general higher than the posterior variances of the gridblocks with measurement information, the multi-scale model transfers information between gridblocks and levels of resolution. As the measurement information
Multi-Scale permeability modeling
11
(a)
10000 0
5000
Iterations
15000
20000
(b)
2
4
6
8
Gridblocks
(c) Fig. 5. Coarsest level of resolution. (a) Version of original field; (b) Posterior mean; (c) Gibbs Sample.
is equally spaced, the transfer of information along the field causes periodic behavior of the posterior variances. Finally, as in the coarsest level of resolution, the posterior variances of the gridblocks on the right end of the field are smaller because of the use of the production data.
7. Conclusions and discussion In this paper, we used a 1-D multi-scale model as a prior for a permeability field. We developed an MCMC algorithm embedded in the Bayesian framework such that a Fluid Flow Simulator can be used at each level of resolution to incorporate the production information. Moreover, the algorithm is performed in a cascade way, from coarser to finer levels. In this way, the realizations at coarser levels guide the realizations at the finer levels leading to a more efficient exploration of the posterior distribution of the multi-scale permeability field. The final result is a multi-scale framework that incorporates data at different resolutions and/or of different types in a coherent and efficient way. We applied the proposed methodology to synthetic data. As expected from intuition, production data gives more information on permeabilities closer to the production well. Finally, the use of data from several different sources allows a very good recovery of the original permeability field. Due to the restriction to 1 dimension, applications to real world problems do not abound.
12
Marco A. R. Ferreira et al.
(a)
10000 0
5000
Iterations
15000
20000
(b)
0
10
20
30
40
50
60
Gridblocks
(c) Fig. 6. Intermediate level of resolution. (a) Version of original field; (b) Posterior mean; (c) Gibbs Sample.
One example of potential application of the proposed methodology is to petroleum fields that exist in the bed of rivers. But the main ideas of the paper survive: the use of a multiscale model defined from coarser to finer levels as prior for the permeability field and the use of the Bayesian paradigm to perform inference about the parameters of interest. Ongoing research is the extension of this work to the 2 dimensional case. The main questions are the development of theoretical multi-scale Markov Random Fields and the development of estimation and simulation procedures. A second step is the adjustment of these theoretical developments to hydrology applications. A. MCMC details The estimation of the parameters of the model is done using an MCMC algorithm. The full conditionals for φ0 and σ02 are trivially obtained as discussed in Chib and Greenberg (1994) and McCulloch and Tsay (1994). The full conditionals for φ1 , σ12 , τ1 , φ2 , σ22 and τ2 are as discussed in Ferreira et al. (2001). Within the MCMC scheme, the various resolution levels of the permeability field are updated with block Metropolis-Hastings steps. In addition, we integrate out the finer levels of resolution when generating a given level of the permeability field. Therefore, the generation of the current level i depends only on the immediate coarser level of resolution
13
0.0
0.002
0.004
0.006
0.008
0.010
Multi-Scale permeability modeling
1
2
3
4
5
6
7
8
0.0015
0.0025
0.0035
0.0045
(a)
0
10
20
30
40
50
60
(b) Fig. 7. Posterior variance across the field. (a) coarse level. (b) intermediate level. Solid line = posterior variance. Dashed line = variance of the measurement data
i − 1. The full conditional for σ21 is equal to IG(n∗ , n∗ s∗ /2), where n∗ = n + T and n∗ s∗ = n s + (pobs − f (X1 ))0 (pobs − f (X1 )). The full conditional for σ20 is obtained analogously. References Chib, S. and Greenberg, E. (1994). Bayes inference in regression models with arma(p,q) errors. Journal of Econometrics, 64:183–206. Craig, P. S., Goldstein, M., Seheult, A. H., and Smith, J. A. (1996). Bayes linear strategies for history matching of hydrocarbon reservoirs. In Bernardo, J. M., Berger, J. O., Dawid, A. P., and Smith, A. F. M., editors, Bayesian Statistics 5, pages 69–95. Clarendon Press, Oxford. Cunha, L. B., Oliver, D. S., Redner, R. A., and Reynolds, A. C. (1996). A hybrid markov chain monte carlo method for generating permeability fields conditioned to multiwell
14
Marco A. R. Ferreira et al.
pressure data and prior information. In SPE Annual Technical Conference and Exhibition, Denver. SPE. Ferreira, M. A. R., West, M., Lee, H., and Higdon, D. (2001). A class of multi-scale time series models. Technical report, Duke University - ISDS. Floris, F. J. T., Bush, M. D., Cuypers, M., Roggero, F., and Syversveen, A.-R. (1999). Comparison of production forecast uncertainty quantification methods — an integrated study. 1st Symposium on Petroleum Geostatistics, Toulouse, 20–23 April 1999. Gamerman, D. (1997). Markov Chain Monte Carlo. Chapman Hall. Gilks, W. R., Richardson, S., and Spiegelhalter, D. J., editors (1996). Markov Chain Monte Carlo in Practice. Chapman Hall. Hegstad, B. K. and Omre, H. (1997). Uncertainty assessment in history matching and forecasting. In E. Y. Baafi and N. A. Schofield, editor, Geostatistics Wollongong ’96, Vol. 1., pages 585–596. Kluwer Academic Publishers. King, M. J. and Datta-Gupta, A. (1998). Streamline simulation: a current perspective. In Situ, 22(1):91–140. Lee, H., Higdon, D., Bi, Z., Ferreira, M. A. R., and West, M. (2002). Markov random field models for high-dimensional parameters in simulations of fluid flow in porous media. Technometrics. To appear. McCulloch, R. E. and Tsay, R. S. (1994). Bayesian analysis of autoregressive time series via the gibbs sampler. Journal of time series analysis, 15:235–250. Neuman, S. P. and Yakowitz, S. (1979). A statistical approach to the problem of aquifer hydrology: 1. theory. Water Resources Research, 15:845–860. Oliver, D. S. (1994). Incorporation of transient pressure data into reservoir characterization. In Situ, 18:243–275. Oliver, D. S., Cunha, L. B., and Reynolds, A. C. (1997). Markov chain monte carlo methods for conditioning a permeability field to pressure data. Mathematical Geology, 29:61–91. Vasco, D. W. and Datta-Gupta, A. (1999). Asymptotic solutions for solute transport: A formalism for tracer tomography. Water Resources Research, 35(1):1–16. Yeh, W. W. (1986). Review of parameter identification in groundwater hydrology: the inverse problem. Water Resources Research, 22:95–108.
