Transcript Document

Trends in quantifying
hydrogeologic uncertainties
Shlomo P. Neuman
Hydrology and Water Resources
University of Arizona
UA Uncertainty Workshop April 25 – 26, 2008
1
Sources of Hydrogeologic Uncertainty
• Incomplete system knowledge (conceptual)
• Measurement/sampling/interpretive (data)
errors
• Unknown/random heterogeneity (properties)
• Unknown/random forcing terms (drivers)
• Disparate data/model scales (conceptual/data)
• The Question: How to account jointly for
these uncertainties?
2
Trends
Use maximum likelihood Bayesian model
averaging (MLBMA) to account jointly for
• Conceptual model uncertainty
• Model parameter uncertainty
• Data errors
Describe multiscale spatial variability
geostatistically
Model flow/transport stochastically
3
Bayesian Model Averaging (BMA)
(Hoeting et al. 1999)
M =  M1 ,..., M K  = set of models considered
 = quantity to predict
K
p  D   p  M k , D p  M k D
k 1
= posterior distribution of  given data D
p   M k , D  = posterior distribution under Mk
p  M k D  = posterior model probability (weight)
All probabilities implicitly conditional on M
4
BMA (2)
Prediction (posterior mean) of :
K
E   D    E   D, M k  p  M k D 
k 1
Predictive uncertainty (posterior variance):
K
Var   D   Var   D, M k  p  M k D 
k 1
within-model variance
K

  E   D, M k   E   D 
k 1
 p  M D
2
between-model variance
k
5
Maximum Likelihood BMA (MLBMA)
(Neuman 2003; Ye et al. 2003 – 2008)
To render BMA computationally feasible set

p   M k , D   p  M k , ˆk , D

 1

exp   KICk  p  M k 
2


p  M k D  K
 1

exp

KIC


l  pMl 
 2

l 1
ˆ
= ML estimate of parameters k
KIC = Kashyap’s model discrimination criterion
k
6
MLBMA (2)
 Obtain ˆk , estimation covariance and Kashyap’s
criterion using ML inversion
• with/without prior information about parameters
• deterministic models (Carrera and Neuman 1986)
• stochastic moment models (Hernandez et al. 2002)
 Use Monte Carlo or stochastic moment model to
estimate predictive uncertainty
E   M k , ˆk , D   E   M k , D 


Var   M k ,ˆk , D   Var   M k , D 


7
Apache Leap Research Site (ALRS)
Unsaturated fractured tuff
]
x [m
40
30
20
20
10
1-m-scale packer tests
y
0
W3
-10
-20
W2 W2A
0
Log10k data
W1
X3
Y3
X2
X1
[m
]
0
Y2
Y1
-20
0
-10
Z1
-10
-20
Z2
-20
-30
Z3
V2
-40
V3
-30
z [m]
z [m]
Geostatistical models
V1
-40
20
y
[m
]
40
30
20
0
10
0
-10
-20
]
x [m
-20
8
Alternative Variogram Models
Fractal
1.8
Power
560
1.6
Homogeneous
1.4
1346
Exponential
Polynomial drift
1st-order
2nd-order
Variogram
Spherical
1.2
2510
1
3338
5312
0.8
4934
4736
0.6
4168
4196
1830
0.4
446
0.2
0
0
5
10
15
20
25
30
Lag Distance (m)
9
Posterior Model Probabilities
Model
Pow0
Exp0
Exp1
Exp2
Sph0
Sph1
Sph2
Parameters
NLL
2
352.1
2
361.0
6
341.6
12
330.4
2
379.1
6
349.6
12
338.9
KIC
369.2
370.1
369.5
416.7
390.5
378.1
424.6
p(Mk)
p(Mk|D)(%)
1/7
35.29
1/7
26.58
1/7
37.61
1/7
0
1/7
0
1/7
0.51
1/7
0
p(Mk|D)(%)
35.29
26.58
37.61
N/A
N/A
0.51
N/A
Sensitive to choice of prior model probabilities
Three models dominate; could discard the rest
No justification for discarding all but one model (as normally done)
10
Z(m)
logk
-14
-14.1
-14.2
-14.3
-14.4
-14.5
-14.6
-14.7
-14.8
-14.9
-15
-15.1
-15.2
-15.3
-15.4
-15.5
-15.6
-15.7
-15.8
-15.9
-16
-16.1
-16.2
-16.3
-16.4
Z(m)
BMA Results (1)
Kriging variance
1.25
1.2
1.15
1.1
1.05
1
0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
(a)
Posterior mean =
Weighted sum of
Kriging estimates
Posterior variance =
Weighted sum of
within-model +
between-model
variances
Weights = Posterior
model probabilities
0
-10
-20
-30
-10
0
10
20
30
40
X(m)
(b)
0
-10
-20
-30
-10
0
10
20
30
40
X(m)
11
BMA Results (2)
1
pow0
exp0
exp1
sph1
BMA
0.9
0.8
Weighted sum of
model probabilities
0.7
Empirical cdf
Posterior probability =
0.6
0.5
0.4
0.3
0.2
0.1
0
-17
-16
-15
-14
-13
log10k
Variogram
model
Variance of 
Pow0
Exp0
Exp1
Sph1
BMA
0.33
0.13
0.47
0.40
0.40
12
Cross-validation
For each of 6 boreholes
Ignore measured log10k data
Estimate variogram parameters and model
probabilities based on remaining data
Assess/compare predictive capabilities of
models & BMA
13
Predictive Log Score
For single model:  ln p( DT | Mk , DB )
K

 
T
B
B

ln
p
(
D
|
D
)


ln
p
D
M
,
D
p
M
D
For BMA:

k
k
T
B
k 1
DT = Predicted values
Pow0
Predictive log score 34.109
DB = Data used
Exp0
Exp1
BMA
35.244 33.968 31.394
Smallest
The smaller the less information is lost
14

Predictive Coverage
Proportion of cross-validation predictions falling
in Monte Carlo simulated 90% prediction interval
Pow0
Exp0
Exp1
BMA
Predictive coverage 86.49
80.83
83.74
87.46
Largest
The larger the better the model’s predictive capability
15
GEOSTATISTICAL
CHARACTERIZATION OF
HIERARCHICAL MEDIA
16
Hierarchical Media
Sedimentary deposits
Ritzi et al. 2006
Fractured rocks
Barton 1995
17
Hierarchical Media Properties are
Multiscale (Nonstationary)
Fractal plots for fracture trace maps (Barton 1995)
18
Permeability of Upper Oceanic Crust
(Fisher 2005)
19
Longitudinal Apparent (Fickian)
Dispersivity (Neuman 1995)
20
Hierarchical Media Properties appear to
fit Traditional Stationary Variograms
Ritzi et al. 2004, Dai et al. 2005, Ritzi & Allen-King 2007 concluded:
Hierarchical sedimentary architecture & ln K tend to have
stationary “exponential-like” transition probabilities & variograms,
the latter arising from the former.
21
How can nonstationary hierarchical
media yield stationary variograms?
Answer: Stationarity is artifact of sampling
over finite domain (window)
Proposed Resolution: Use truncated power
variograms (Di Federico & Neuman
1997; Di Federico et al. 1999) or TPVs
TPV = stationary variogram of nonstationary
fractal field sampled on given support
scale over finite window
22
Fractal Decomposition
Power variogram
  s   Co s
2H
is weighted integral
of exponential or
Gaussian modes:

dn
  s      s, n 
n
0
n  1/   mode number
  integral scale
23
Fractal Decomposition
Introducing upper and lower cutoffs:
nu
dn
  s; nl , nu      s; n 
n
nl
Lower cutoff renders field statistically
homogeneous.
Sample variograms fitting standard models
of homogeneous fields often fit   s; nl , nu 
equally well.
24
1
0.5
0
0
2
4
Normalized lag
H = 0.375
1
0.5
0
0
Normalized variograms
H = 0.25
Normalized variograms
Normalized variograms
TPVs are hard to distinguish from
traditional stationary models
2
4
Normalized lag
H = 0.55
1
0.5
0
0
2
4
Normalized lag
Exponential (solid) compared with G-TPV (dashed)
G-TPV = TPV based on Gaussian modes
H = Hurst scaling coefficient
25
Tübingen Case Example
5 pumping tests in B wells
interpreted graphically
(Neuman et al. 2004, 2007):
2
TG  2.2  10 m / s
2
Integral Scale (ln T) = 2.5 m
Variance (ln T) = 1.5
Support = 0 Window = 25 m Note IS/W = 1/10!
26
Tübingen Case Example
312 flowmeter K values from
all B and F wells were locally
upscaled to obtain
2
TG  2.4  10 m / s
2
Variance (ln T) = 1.4
Window = 306 m
Insufficient data to estimate integral scale
27
Tübingen Case Example
TPV allows interpreting
these two sets of results
jointly:
For local support = 1 m
obtain
Integral Scale = 50 m
H = 0.08
Local supports of 0.7 and 2 m
yield very similar results
28
Tübingen Case Example
TPV allows cokriging multiscale data.
Using 1-m scale flow meter data to predict
25-m scale Y = ln T (not via block kriging!):
B3
B1B2
300
300
B4
F0
lnT
-2.2
-2.3
-2.4
-2.5
-2.6
-2.7
-2.8
-2.9
-3
-3.1
-3.2
-3.3
-3.4
-3.5
-3.6
-3.7
-3.8
-3.9
-4
-4.1
-4.2
-4.3
-4.4
-4.5
-4.6
250
F1
200
y [m]
F3
150
F6
100
F5
50
F2
B5
2
 lnT
B4
F0
3.8
3.6
250
3.4
F1
3.2
3
2.8
200
2.6
F3
y [m]
F2
B3
B1
B2
B5
2.4
2.2
2
150
1.8
1.6
1.4
1.2
F6
100
1
0.8
0.6
F5
50
F4
50
100
150
F4
200
250
x [m]
Y estimates
300
50
100
150
200
250
300
x [m]
Estimation variance 29
Tübingen Case Example
Using 1-m scale flow meter data and B-well
test results jointly to predict 1-m scale
Y = ln T by cokriging:
B3
B1
300
F2
B3
B5
B2
F0
B1
300
F2
B4
B5
B2
F0
B4
-1
250
200
-2
F3
-2.5
-3
150
1.4
y [m]
y [m]
F1
1.6
250
-1.5
F1
1.2
200
F3
1
0.8
150
-3.5
F6
-4
100
0.6
F6
100
0.4
-4.5
F5
-5
50
F5
0.2
50
0
F4
50
100
150
F4
200
250
x [m]
Y estimates
300
50
100
150
200
250
300
x [m]
Estimation variance30
Anomalous Transport in
Randomly Heterogeneous Media
31
Contaminant Plumes in
Groundwater at Hanford (2005)
32
Classic Fickian Model
for Passive Tracer
• Dispersive flux =
c
  d m   L,T v 
 xL,T
d L ,T
d m ,  L ,T
= constants
• Advection-dispersion equation (ADE)
c
   dc  vc   g 
t
33
Classic Fickian Model
Plume is Gaussian
Its variance grows linearly with time
34
Heterogeneity
(on all scales)
tends to render
dispersion
non-Fickian
(anomalous)
Photo by John Selker, Oregon State Univ.,
from a trench face near the 200 East Area,
Hanford Site, Washington
35
Field manifestation of
non-Fickian dispersion
At Borden dispersivity of solute slug
under mean-uniform flow varies:
• Linear in early travel-time/meandistance (1:1 log-log slope)
• Longitudinal stabilizes at “Fickian”
asymptote
• Transverse peaks, then diminishes
toward (Zhang and Neuman 1990;
Dentz et al. 2002; Attinger et al.
2004)
– Local value or 0 in 2D
– Larger nonzero asymptote in 3D
(Zhang and Neuman 1990)
Localized stnL models
36
Laboratory manifestation of
non-Fickian dispersion
Similar behavior is exhibited in 3D matched index
particle tracking velocimetry experiments of Moroni et al.
(2007) in terms of normalized longitudinal (left) and
transverse (right) dilution indexes:
37
Simulated (localized stnADE) Effects of
Boundaries & Conditioning on D
(Morales-Casique et al. 2005)
38
Field + lab manifestations of
non-Fickian dispersion
Preasymptotic scaling also
observed in radial field test
(Peaudecerf and Sauty 1978)
and 2D lab tests (Silliman
and Simpson 1987; Sternberg
et al. 1996)
Apparent dispersivities
worldwide increase with
scale at supralinear rate
(Neuman 1990, 1995)
excluding calibrated models
(diamonds)
39
Laboratory manifestations of
non-Fickian dispersion
Anomalous early and late (tailing)
BTCs in Berea Sandstone core
(right) and Aiken clay loam column
(left); CTRW solutions fitted by
Cortis and Berkowitz (2004)
40
Pore-scale simulation of
non-Fickian dispersion
(Zhang & Lv 2007)
41
Space-time nonlocal theory based
on local ADE (stnADE)
(Neuman 93 - Morales et al. 05)
Premises:
• All variables defined on support scale 
• K(x) is correlated random field
• Therefore v(x,t) = random function
c
• Hence ADE
   dc  vc   g 
t
d =  - scale dispersion tensor
g = random source
is stochastic.
42
Monte Carlo Solution
•
•
•
•
Requires high-resolution grid
Yields multiple random images of c
Realistic but equally likely (nonunique)
Averaging them yields (at the least)
C = (conditional) ensemble mean of c
= optimum unbiased predictor of c
= deterministic = smooth relative to c
= (conditional) ensemble variance of c
 c2 = measure of predictive uncertainty
43
stnADE Transport Equation
C
   dC  Q  VC   G 
t
V = (conditional) mean v
Q = (conditional) mean dispersive flux
C, V, G = averages over MC realizations
C, V, G, Q defined on scale 
44
stnADE Transport Equation (2)
• Q = integro-differential (nonlocal)
Q( x, t )     ( x, t, y, s)C ( y, s) ds dy + source / bdry terms
y s
 = (conditional) nonlocal parameter
• stnADE provides second moments (predictive
covariance) of C and solute flux
45
Computational Example
46
Transport Conditions

c
0
x2
Boundary conditions
c
v1c  Dd
0
x1
c
0
x1
x2
x1

Initial concentration
c
0
x2
47
MC (left) and moment equation (right) velocity statistics
48
2nd order (left) and iterative 2+ order (right) C
49
Longitudinal Mean Mass Flux
50
Lagrangian Motions (stnL)
(Cushman & Ginn 1993 – Cushman
1997)
• Solute particles move by
– Brownian motion (local dispersion, left out below)
– Correlated random velocity v(t) = dX(t)/dt
X(t) = random particle displacement
• v(t) or increments stationary
• Dispersion flux of particle probability p(x,t):
t
Q( x, t )  
 D ( y, t, )
2
x y
p( x - y, t   ) dy d
0 
• stnL = stationary (convolution) stnADE version
51
Fractional ADE (fADE)
(Benson 1998 – Zhang et al. 2007)
• No underlying process assumptions
• Deterministic fADE postulated directly
• Most common fADE is 1D:

C
C
 C
1   2
 V
D 
t
x
x
• Generalized to multiple dimensions,
fractional temporal terms, spatially varying
V(x) and D(x)
52
fADE with variable V and D
Integrodifferential equivalents contain
dispersive flux (convolution)
t 
C  x  y , t   
Q  x, t       x, t; y , 
dy d
x  y
0 
with (time-localized!) kernels
D  x     H  y 
D  x  y     H  y 


 1
 2    y
  2    y  1
FF-ADE
FD-ADE, FFD-ADE
53
fADE with variable V and D (3)
• Dirac delta renders Q nonlocal in space
but local in time
• fADE with variable V and D is therefore
not equivalent to stnADE (as proposed by
Zhang et al. 2007)
54
Continuous Time Random Walk
(CTRW) (Berkowitz et al. 2006)
• Solute particles at given concentrations
hop between nodes-sites of fixed grid
• Complete instantaneous mixing of solute
at each node-site
• Stochastic master (balance) equation:
c  x, t 
   w  y , x  c  x, t    w  y , x  c  y , t 
t
y
y
w(y,x) = transition rate from y to x
55
Continuous Time Random Walk
(CTRW) (2)
• If w(y,x) = incoherent stationary random
field then generalized master equation:
t
C  x, t 
    y  x, t    C  x, t  d
t
y 0
t
    x  y, t    C  y, t  d
y
0
  z,  
  z,   
1  s   
 s       z,  
z
56
Continuous Time Random Walk
(CTRW) (3)
  z, s 
  z,   = Laplace transform of
  z, s  = probability rate of displacements
z over time intervals s
 s  s  = marginal probability rate of all
such displacements
• Incoherence = lack of statistical
interdependence
• Highly limiting (no spatial correlations)
• Overlooked in hydro CTRW literature
57
Continuous Time Random Walk
(CTRW) (4)
• Coherence would require  z, s  to be
multivariate
• In limit of small displacements obtain pde
for C with dispersive flux
t
Q  x, t    Dl  t   C  x, t  d
0
• Same pde / temporal - convolution obtained
from stnL / stationary stnADE upon
ignoring spatial dependencies!
58