Transcript Stochastic inverse modeling under realistic prior model constraints with multiple-point geostatistics Jef Caers Petroleum Engineering Department Stanford Center for Reservoir Forecasting Stanford, California, USA.

Stochastic inverse modeling under
realistic prior model constraints with
multiple-point geostatistics
Jef Caers
Petroleum Engineering Department
Stanford Center for Reservoir Forecasting
Stanford, California, USA
Acknowledgements
I would like to acknowledge the contributions
of the SCRF team, in particular Andre Journel
and all graduate students
who contributed to this presentation
Quote
“Theory should be as simple as possible,
but not simpler as possible….”
Albert EINSTEIN
Overview
• Multiple-point geostatistics
• Why do we need it ?
• How does it work ?
• How do we define prior models with it ?
• Data integration
• Integration of multiple types/scales of data
• Improvement on traditional Bayesian methods
• Solving general inverse problems
• Using prior models from mp geostatistics
• Application to history matching
Part I
Multiple-point Geostatistics
Limitations of traditional geostatistics
1
2
3
Variograms EW
Variograms NS
1
1.2
1.2
2
0.8
3
0.4
0.8
0.4
0
10
20
30
40
0
10
20
30
40
2-point correlation is not enough to characterize connectivity
A prior geological interpretation is required
and it is NOT multi-Gaussian
Stochastic sequential simulation
• Define a multi-variate (Gaussian) distribution over the
random function Z(u)
• Decompose the distribution as follows
Pr(Z(u1 )  z1,
, Z(uN )  z N ) 
Pr(Z(u1 )  z1 ). Pr(Z(u2 )  z 2 | z1 ). Pr(Z(u3 )  z 3 | z 2 ,z1 )
Pr(Z( uN )  z N | z N 1 ,
,z1 )
Or in its conditional form
Pr(Z(u1 )  z1,
, Z(uN )  z N | (n)) 
Pr(Z(u1 )  z1 | (n)). Pr(Z(u2 )  z 2 | z1,(n)). Pr(Z( u3 )  z 3 | z 2 ,z1, (n))
Pr(Z(uN )  z N | z N 1,
,z1,(n))
Practice of sequential simulation
B1
A
B2
B3
B={B1,B2,B3}
P(A|B) = N(m,s)
m,s given by kriging, depend on
autocorrelation (variogram) function
Multiple-point Geostatistics
Reservoir
= well data
multiple-point
data event
Sequential
simulation
P(A|B)?
A
B
Extended Normal Equations
Attribute S taking K possible states {s k , k  1,
, K}
At location u, unknown binary random variable
1 if S( u)  s k
Ak  
if not
0
u
Pr(A k  1 | S( u  )  s  , ) ?
1 if S( u  )  s  ,   1,
D
if not
0
u
,n
n
  A
1
Traditional kriging : 2 point statistics
n
Pr(A k  1 | S( u  )  s  , )  E[A k ]     (1  E[A  ])
1
Extensions to three point statistics
n
Pr(A k  1 | S( u  )  s  , )  E[A k ]     (1  E[A  ])    (1  E[A A  ])
1

Single Normal Equation
In case of (n+1) - point statistics
Pr(A k  1 | S( u  )  s  , )  Pr(A k  1 | D  1)  E[A k ]  (1  E[D])
.Var[D]  Cov[A k , D]

E[A k D]  E[A k ]E[D]
E[D](1  E[D])
such that
E[A k D]  E[A k ]E[D]
E[D](1  E[D])
E[A k D]  E[A k ]E[D]
Pr(A k  1 | D  1)  E[A k ] 
E[D]
E[A k D] Pr[A k  1, D  1]


E[D]
Pr[D  1]
Pr(A k  1 | D  1)  E[A k ] 
Bayes Rule !
The training image module
Training image module =
standardized analog model
quantifying geo-patterns
P(A|B)=1/4
SNESIM algorithm
A = Mud
Recognizing P(A|B) for all possible A,B
The SNESIM algorithm
Training image
Data template
(data search neighborhood)
Search tree
Construction requires scanning
training image one single time
14 11
5 7
5 3
3 1
1 0
2 5
1 1
1 1
0 3
0 2
3 0
1 0
1 0
1 0
1 0
1 1
2 0
2 0
1 1
0 1
Minimizes memory demand
Allows retrieving all training cpdf’s
for the template adopted!
Probabilities from a Search Tree
Training image
Search neighborhood
5
4
3
2
j=1
1
4u2
3
Search tree
i =1 2 3 4 5
Level 0 (no CD)
u
1
u
Level 1 (1 CD)
.
.
.
.
.
.
1
u2
1
3
0
Level 4 (4 CD)
5
1
u2
1
1
1
1
u
7
1
u2
1
0
1
u
14 11
1
u
2
5
3
0 2
1
u2
1
1
0
1
u2
0
1
1
0
3
2
2
0
1
u
3
1
u2
0
0
5
1
u2
0
1
0
1
1
1
1
1
u2
1
1
u 2 u2
3
3
1
4u2
3
Example
True image
400 sample data
CPU
Training image
Realization
2 facies,
1 million cells
= 4’ 30”
On 1GHz PC
Where do we get a 3D TI ?
2- D Reference Data
500. 000
North
sand
shale
0. 0
0. 0
East
500. 000
Valid training image
Not Valid
Training image requires "stationarity"
Only patterns = "repeated multipoint statistics" can be reproduced
Modular training image
Training Image
Models generated with snesim
using the SAME training image
2- D Reference Data
500. 000
North
sand
shale
0. 0
0. 0
East
500. 000
Modular ?
* no units
* rotation-invariant
* affinity-invariant
Properties of training image
Required
• Stationarity: patterns by definition repeat
• Ergodicity: to reproduce long range feature => large image
• Limited to 4-5 categories
Not required
• Univariate statistics need not be the same as actual field
• No conditioning to ANY data
• Affinity/rotation need not be the same
Part II
Multiple-point Geostatistics
and data integration
Simple question, difficult problem…
The essential data integration problem…
A geologist believes based on geological data that there is 80%
chance of having a channel at location X
P(A|B)
A geophysicist believes based on geophysical data that there is
75% chance of having a channel at location X P(A|C)
A petroleum engineer believes based on engineering data that
there is 85% chance of having a channel at location X
P(A|D)
What is the probability of having a channel at X ?
P(A|B,C,D)?
Combining sources of information
P(A | B,C)  (P(A), P(A | B), P(A | C))  (P(A), P(B | A), P(C | A))
Desirable properties of , 
1. P(A | B,C)  [0,1]
2. Closure P(A | ( ))  P(A | ( ))  1
3. C non-informative of A : P(A | C,( ))  P(A | ( ))
4. C fully informative of A : P(A | C,( ))  1
5. C fully informative of A not occurring : P(A | C,( ))  1
Conditional independence
P(B,C | A)  P(B | A)  P(C | A))
P(A)  P(B | A)  P(C | A)
 P(A | B,C) 
P(B,C)
Requires McMC to eliminate P(B,C)
1. P(A | B,C)  [0,1] ?
O = In practice not necessarily
2. Closure P(A | ( ))  P(A | ( ))  1 ?
3. C non-informative of A : P(A | C,( ))  P(A | ( )) ?
4. C fully informative of A : P(A | C,( ))  0 ?
5. C fully informative of A not occurring : P(A | C,( ))  1 ? YES
Correcting conditional independence
2. Closure P(A | B,C)  P(A | B,C)  1 
P(A)  P(B | A)  P(C | A)  P(A)  P(B | A)  P(C | A)  P(B,C)
Solution : Standardize the expression
P(A)  P(B | A)  P(C | A)
P(A | B,C) 
P(A)  P(B | A)  P(C | A)  P(A)  P(B | A)  P(C | A)
P(A | B)  P(A | C)
P(A)

P(A | B)  P(A | C) P(A | B)  P(A | C)

P(A)
P(A)
This expression honors all conditions 1 - 5
Permanence of ratios hypothesis
P(A | B)  P(A | C)
P(A)
P(A | B,C) 
P(A | B)  P(A | C) P(A | B)  P(A | C)

P(A)
P(A)
can be reshaped as follows
1
a
x b ca
P(A | B,C) 

or

1  x a  bc
b
a
with
1  P(A)
1  P(A | B)
1  P(A | C)
1  P(A | B,C)
a
b
c
x
P(A)
P(A | B)
P(A | C)
P(A | B,C)
xb ca

is a permanence of ratios hypothesis
b
a
Advantages of using ratios
• No term P(B,C), hence McMC is not required
• Work with P(A|B),P(A|C), more intuitive than P(B|A),P(C|A)
• Verifies all consistency conditions by definition
• It is still a form of independence,
Yet dependence can be reintroduced
x b ca
x c

 
b
a
b a
reintroducing

dependence
x c
 
b a
(B,C)
Simple problem…
P(A|B) = 0.80
1
a2
P(A | B,C, D) 
 2
1  x a  bcd
P(A|C) = 0.75 => P(A|BC)
Suppose P(A) = 0.5 => P(A|BC) = 0.92
Suppose P(A) = 0.3 => P(A|BC) = 0.95
= compounding of events
Lesson learned : if geologist and geophysicist agree
for almost 80%, you can be even more certain that there
is a channel !
Example reservoir
Training image
P(A|B)
Realization
P(A|C)
ence2-D
Data
Reference Data
226.000
.000
1.000
SingleRealization
realization
226.000
1.000
North
North
North
shale
sand
0.0
0.0
ast
0.0
0.0
0.0
0.0 112.000
East
112.000
0.0
East
112.000
0.0
East
112.000
P(A|C), A = single-point !
Realization
Realization
P(A|C)
2-D Reference Data
226.000
.000
1.000
North
North
?
When combing P(A|B) from
geology and1.000
P(A|C) from seismic
to P(A|BC), 0.900
0.8000
‘A’ is still a0.6000
single point event !
0.7000
0.5000
0.4000
 Certain patterns,
such as
local rotation0.3000will be ignored
0.2000
0.1000
0.0
0.0
0.0
0.0
East
112.000
0.0
East
112.000
 Honor seismic
only as a
0.0
single-point probability
Concept of MODULAR training image
Modular Training Image
2- D Reference Data
500. 000
North
sand
shale
0. 0
0. 0
East
500. 000
Modular ?
* Stationary patterns
* rotation-invariant
* affinity-invariant
* no units
Models generated with snesim
using the SAME training image
Local rotation angle from seismic
P(A|C)
Local angle
2-D Reference Data
226.000
220
1.000
200
180
160
North
140
120
100
80
60
40
0.0
20
0.0
0.0
East
112.000
0
20
40
60
80
100
Results
2 Realization
realizations with
angle
Realization
226.000
without
angle
Realization
226.000
226.000
220
1.000
200
0.900
180
0.8000
shale160
0.7000
0.6000
120
0.5000
100
0.4000
80
0.3000
60
0.2000
North
North
North
140
sand
0.1000
40
0.0
20
0.0
0.0
0.0
0.0
East
112.000
0.0
East
112.000
0
20
40
60
80
100
0.0
East
112.000
Constrain to local “channel features”
2-DP(A|C)
Reference Data
Feature map
226.000
226.000
Realization
226.000
1.000
sand
North
North
North
shale
not determin
0.0
0.0
0.0
0.0
East
112.000
Hard data
from seismic
0.0
East
112.000
Soft data
from seismic
0.0
0.0
East
112.000
Part III
Inverse modeling with
multiple-point geostatistics
Application to history matching
Production data does not inform
geological heterogeneity
Geologist 1
A Petroleum Engineer
Realization matched to Fractional Flow
Realization matched only to Fractional Flow Data
50.0
1000.0
800.0
800.0
600.0
600.0
400.0
North
North
50.0
1000.0
400.0
1
1.0
a
a
a
a
a
a
aa
a a
aa
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
Initia l P erm e ab
a
200.0
a
a
a
200.0
a
a
a
a
a
a
a
aa
a
a
a
aa
0.0
0.0
East
0.0
50.0
?
Disagreeing Geologist
Realization matched to Fractional Flow
50.0
0.0
East
50.0
Geologist 2
The Reference Permeability Field
50.0
1000.0
0.8
Water Cut
0.0
0.0
Water cut
a
a
1000.0
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
0.6
a
a
a
a
a
a
a
a
a
0.4
a
a
a
a
0.2
a
800.0
800.0
600.0
600.0
400.0
North
North
a
400.0
0.0
a
a
a
a
a
a
a
0
0
0.0
East
50.0
0.0
East
50.0
600
800
1000
Time (days) 1000
a
0.0
0.0
400
tim e , da ys
0.0
0.0
a
a
200
200.0
200.0
a
Approach
Methodology
Define a non-stationary Markov chain that moves a
realization to match data, two properties
• At each perturbation we maintain geological realism
use term P(A|B)
• Construct a soft data set “P(A|D)” such that we move
the current realization as fast as possible to match the
data => Optimization of the Markov chain at each step
Methodology: two facies
Some notation:
1

I(u)  
0

if facies s1 occurs
u  (x, y,z)
if facies s0 occurs
D = set of historic production data (pressures, flows)
Initial guess realization:
i(o) (u) u
Realization at iteration ( )
i (u) u
( )
Define a Markov chain
i( ) ( u) u
?
i(
1)
( u) u
Define a transition matrix:
Pr{I (
1)
( u)  1 | D,i ( ) ( u)  0}  ?
Pr{I (
1)
( u)  0 | D,i ( ) ( u)  1}  ?
Pr{I (
1)
( u)  1 | D,i ( ) ( u)  1}  ?
Pr{I (
1)
( u)  0 | D, I ( ) ( u)  0}  ?
Transition matrix
2 x 2 transition matrix describes the probability of changing
facies at location u and we define it as follows
Pr{I (
1)
( u)  1 | D,i ( ) ( u)  0}  rD Pr{I( u)  1}
Pr{I (
1)
( u)  0 | D,i ( ) ( u)  1}  rD (1  Pr{I( u)  1})
Pr{I (
1)
( u)  1 | D,i( ) ( u)  1}  1  rD (1  Pr{I( u)  1})
Pr{I (
1)
( u)  0 | D,i ( ) ( u)  0}  1  rD Pr{I( u)  1}
rD  [0,1]
rD is a constant over the reservoir and depends on production D
Parameter rD
A  {I (
1)
1)
Pr{I (
( u)  1}
( u)  1 | D,i ( ) ( u)  0}  rD Pr{I( u)  1}
rD  [0,1]

P(A | D) = rDP(A) if i (l) (u) = 0
1)
Pr{I (
( u)  1 | D,i ( ) ( u)  1}  1  rD Pr{I( u)  0}
rD  [0,1]

P(A | D) = 1 - rD (1 - P(A)) if i (l) (u) = 1
 P(A | D) = (1- rD ) ×i (l) (u) + rD ×P(A)
Determine rD
Use P(A|D) as a probability model
in multiple-point geostatistics
 Combine P(A|B) (from training image) with P(A|D)
from production data D into P(A|B,D)
 Allows generating iterations that are consistent with prior
geological vision
 Allows combining geological information with production
data
 Allows determining an optimal value for rD as follows…
rD determines a “perturbation”
rD=0.01
Some initial model
initial model
training image
r_D = 0.05
50.000
50.000
150.000
facies 1
North
North
facies 1
facies 0
0.0
facies 0
0.0
East
50.000
0.0
rD=0.1
East
50.000
rD=0.2
r_D = 0.1
North
0.0
r_D = 0.2
50.000
50.000
facies 1
North
North
facies 1
facies 0
0.0
facies 0
0.0
0.0
East
50.000
0.0
East
50.000
0.0
rD=0.5
0.0
rD=1
r_D = 0.5
50.000
facies 1
North
North
facies 1
facies 0
0.0
facies 0
0.0
East
150.000
r_D = 1
50.000
0.0
East
50.000
0.0
East
50.000
Find rD that matches best
the production data
= one-dimensional optimization
Complete algorithm
• Construct a Training Image with the desired geological
continuity constraint
• Use snesim (P(A|B)) to generate an initial guess i
(o)
(u)
• Until adequate match to production data D
• Define a soft data P(A|D) as function of rD
• Perform snesim with P(A|D) to generate a new guess
i(rD1) (u)
• Find the value of rD that matches best the data D
training image
150.000
North
Examples
0.0
0.0
Reference model
50.000
North
facies 0
0.0
0.0
East
150.000
P
facies 1
I
East
50.000
Generate 10 reservoir models
that
1. Honor the two hard data
2. Honor fractional flow
3. Have geological continuity
similar as TI
0.45
reference
match
init
0.4
Single model
0.35
0.3
Initial model
Iteration 1, r_D=0.21
50.000
0.25
fw
50.000
0.2
facies 1
facies 1
North
North
0.15
0.1
facies 0
facies 0
0.05
0.0
0
0
0.0
0.0
East
50.000
0.0
East
5
10
15
20
25
timesteps (days)
50.000
30
35
40
1
0.9
Iteration 5, r_D=0.50
50.000
50.000
facies 1
North
North
facies 1
Objective function
0.8
Iteration 3, r_D=0.52
0.7
0.6
0.5
0.4
0.3
facies 0
facies 0
0.2
0.1
0.0
0.0
0.0
East
50.000
0.0
East
50.000
0
0
1
2
3
4
5
Iterations
6
7
8
9
0.65
0.6
Iteration 7, r_D=0.31
0.55
Iteration 9, r_D=0.24
50.000
50.000
0.5
value of r
facies 1
North
North
facies 1
D
0.45
0.4
0.35
facies 0
facies 0
0.3
0.25
0.0
0.2
0.0
0.0
East
50.000
0.0
East
50.000
0.15
1
2
3
4
5
Iterations
6
7
8
9
rD values,
single 1D optimization
It er at i on1
It er at io n3
It er at i on 2
0.03
0.024
0.022
0.02
Objective function
0.022
0.025
0.02
0.018
0.02
0.016
Objective function
0.015
0.014
0.012
0.01
0.01
0.005
0
0.5
1
0.008
0.5
- 3
It er at i on4
10
0.018
9
0.016
8
0.014
7
0.012
6
0.01
5
0.016
0.014
0.012
0.01
0.008
0
0.02
0.018
x 10
1
0.006
0
It er at i on 5
0.5
1
It er at io n6
0.026
0.024
0.022
0.02
0.018
0.016
0.014
0.008
0
0.5
- 3
5
x 10
1
12
11
4
10
3.5
9
3
8
2.5
7
2
6
1.5
5
0
0.5
t - valu e
0.5
- 3
It er at i on7
4.5
1
0.012
4
0
1
x 10
4
0
1
0.01
0
It er at i on 8
0.5
1
It er at io n9
0.012
0.01
0.008
0.006
0.004
0.002
0.5
t- v al ue
rD value
1
0
0
0.5
t - v al ue
1
North
North
fac ies 0
0.0
facies 0
Different geology
0.0
0.0
50.000
East
0.0
East
50.000
0.5
50.000
50.000
50.000
0.45
Iteration
Iteration7 1
50.000
50.000
0.900
0.900
0.8000
fac ies 1
0.8000
facies 1
0.7000
0.7000
0.6000
0.6000
0.5000
0.4000
facies 0
1.000
North
North
North
North
North
facies 1
reference
match
init
0.4
1.000
0.35
0.3
fw
Iteration
6
Initial
model
Reference
model
0.25
0.2
0.5000
0.15
0.4000
0.3000
fac ies 0
0.3000
facies 0
0.2000
0.2000
0.1
0.05
0.1000
0.1000
0.0
0.0
0.00.0
0.0
0.0
East
East
50.000
East
East
Iteration 3
1.000
1.000
0.900
0.900
0.900
0.8000 0.8000
0.8000
0.7000 0.7000
0.7000
0.6000 0.6000
0.6000
0.4000
0.5000
0.4000
North
North
1.000
0.5000
0.2000
0.1000
0.0
0.0
0.0
0.0
East
East
10
15
20
25
timesteps (days)
30
35
40
reference
match
init
50.000
0.5000
0.4000
0.3000 0.3000
0.0
5
50.000
50.000
0.35
Reference
model 2
Iteration
50.000
North
50.000
East
50.000
50.000
50.000
0.3
0.25
0.2
fw
0.0
0
0
0.0
0.0
0.0
0.0
0.0
0.15
0.3000
0.2000
0.2000
0.1000
0.1000
0.0
0.0
0.1
0.05
0.0
50.000
0.0
East
50.000
0
0
5
10
15
20
25
timesteps (days)
30
35
40
Initial model
100.000
More wells
Iteration 6
fac
No rth
100.000
facies 1
facies 1
Reference model
100.000
No rth
fac
facies 1
facies 0
0.0
0.0
North
facies 0
East
Well 1
Well 2
0.8
0.0
0.0
East
0.7
100.000
facies 0
100.000
0.35
reference
match
init
0.3
0.6
reference
match
init
0.25
0.5
East
100.000
0.15
100.000
0.3
0.0
0.0
fw
fw
0.2
0.4
0.1
0.2
0.05
0.1
0
0
100.000
10
20
30
timesteps (days)
Well 3
facies 1
0.5
facies 1
0.4
0
0
40
No rth
Iteration 8
Iteration 4
fac
10
20
30
timesteps (days)
40
Well 4
0.7
reference
match
init
0.6
reference
match
init
fac
No rth
0.5
0.3
facies 0
facies 0
fw
fw
0.4
0.2
0.3
0.0
0.0
0.2
East
100.000
0.1
0.1
0
0
0.0
0.0
East
100.000
10
20
30
timesteps (days)
40
0
0
10
20
30
timesteps (days)
40
Hierarchical matching
* First choose fixed permeability per facies,
perturb facies model
* Then, for a fixed facies
perturb the permeability within facies
(using traditional methods, ssc, gradual deformation
Example
Reference model
50.000
P
North
facies 1
facies 0
I
0.0
0.0
East
50.000
Initial guess
Iteration 1
50.000
50.000
K=500mD
No rth
No rth
K=500mD
K=50mD
0.0
0.0
K=50mD
0.0
Klow = 50
Khigh = 500
50.000
East
0.0
Klow = 150
Khigh = 750
Klow = 50
Khigh = 500
East
50.000
Reference model
Iteration 2
50.000
Iteration 3
50.000
50.000
facies 1
North
K=729mD
No rth
No rth
K=729mD
K=12mD
0.0
0.0
Klow = 12
Khigh = 729
0.0
50.000
East
facies 0
K=12mD
0.0
Klow = 12
Khigh = 729
East
50.000
0.0
0.0
East
Iteration 4
50.000
No rth
K=694mD
Klow = 11
Khigh = 694
plot5.ps
K=11.5mD
0.0
0.0
East
50.000
Tue Oct 23 14:13:10 PDT 2001
Results
50.000
Results
Initial Guess, k_1=500, k_0=50
50.000
North
0.4
fw
0.3
0.0
0.0
East
0.4
reference
match
init
0.3
0.2
0.2
0.1
0.1
reference
match
init
50.000
After iter 1, k_1=721, k_0=12
50.000
10
20
30
0
40
10
20
30
40
30
40
North
0
0.4
0.3
0.4
reference
match
init
0.3
reference
match
init
0.0
fw
East
50.000
fw
0.0
0.2
0.2
After iter 4, k_1=694, k_0=11.5
50.000
0.1
North
0
0.1
10
20
30
timesteps (days)
0.0
0.0
East
50.000
40
0
10
20
timesteps (days)
More realistic
Reference
matched model
Initial model
2-D Reference Data
2-D Reference Data
2-D Reference Data
50.000
50.000
0.0
400.000
400.000
400.000
300.000
300.000
300.000
200.000
200.000
100.000
100.000
100.000
0.0
0.0
0.0
North
North
North
50.000
200.000
0.0
0.0
0.0
East
50.000
50.000
0.0
East
0.0
50.000
North
facies 0
East
50.000
facies 1
facies 1
North
North
facies 1
0.0
50.000
50.000
50.000
0.0
East
facies 0
facies 0
0.0
0.0
0.0
East
50.000
0.0
East
50.000
Conclusions
What can multiple-point statistics provide
• Large flexibility of prior models, no need for math. def.
• A fast, robust sampling of the prior
• A more realistic data integration approach than traditional
Bayesian methods
• A generic inverse solution method that honors prior
information
More on conditional independence
P(B | A) ?
B  f B (A)   B similarly C  f C (A)  C
Conditional independence =  B , C are independent
P(B  b | A,C  c)  P(  B  b  f B (A) | C  f C (A)  c, A)

P(  B  b  f B (A) | A), b,c

P(B | A)
Q? why should B and C be independent unless they are
homoscedastic, i.e. independent of A ?
Q? Is it not a mere transfer of independence hypothesis to
B and C