Trends in quantifying hydrogeologic uncertainties

1. Trends in quantifying hydrogeologic uncertainties Shlomo P. Neuman Hydrology and Water Resources University of Arizona UA Uncertainty Workshop April 25 – 26, 2008

2. Sources of Hydrogeologic Uncertainty • Incomplete system knowledge (conceptual) • Measurement/sampling/interpretive (data) errors • Unknown/randomheterogeneity (properties) • Unknown/random forcing terms (drivers) • Disparate data/model scales (conceptual/data) • The Question:How to account jointlyforthese uncertainties?

3. Trends • Use maximum likelihood Bayesian model averaging (MLBMA) to account jointly for • Conceptual model uncertainty • Model parameteruncertainty • Data errors • Describe multiscale spatial variability geostatistically • Model flow/transport stochastically

4. Bayesian Model Averaging (BMA)(Hoeting et al. 1999) M== set of models considered • = quantity to predict = posterior distributionof  given dataD = posterior distribution under Mk = posterior model probability (weight) All probabilities implicitlyconditional onM

5. BMA (2) Prediction (posterior mean)of: Predictive uncertainty (posteriorvariance):

6. Maximum LikelihoodBMA (MLBMA)(Neuman 2003; Ye et al. 2003 – 2008) To render BMA computationally feasible set = ML estimate of parameters k KIC = Kashyap’s model discrimination criterion

7. MLBMA (2) • Obtain , 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) • UseMonte Carloorstochastic moment model to estimate predictive uncertainty

8. Apache Leap Research Site (ALRS) • Unsaturated fractured tuff • 1-m-scale packer tests • Log10k data • Geostatistical models

9. Alternative Variogram Models • Fractal • Power • Homogeneous • Exponential • Spherical • Polynomial drift • 1st-order • 2nd-order

10. Model Pow0 Exp0 Exp1 Exp2 Sph0 Sph1 Sph2 Parameters 2 2 6 12 2 6 12 NLL 352.1 361.0 341.6 330.4 379.1 349.6 338.9 KIC 369.2 370.1 369.5 416.7 390.5 378.1 424.6 p(Mk) 1/7 1/7 1/7 1/7 1/7 1/7 1/7 p(Mk|D)(%) 35.29 26.58 37.61 0 0 0.51 0 p(Mk|D)(%) 35.29 26.58 37.61 N/A N/A 0.51 N/A Posterior Model Probabilities • 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)

11. BMA Results (1) • Posterior mean = Weighted sum of Kriging estimates • Posterior variance = Weighted sum of within-model + between-model variances • Weights = Posterior model probabilities

12. Variogram model Pow0 Exp0 Exp1 Sph1 BMA Variance of  0.33 0.13 0.47 0.40 0.40 BMA Results (2) Posteriorprobability = Weighted sumofmodelprobabilities

13. Cross-validation • For each of 6 boreholes • Ignore measured log10k data • Estimate variogram parametersand model probabilities based on remaining data • Assess/compare predictive capabilities of models & BMA

14. Predictive Log Score Forsingle model: ForBMA: DT =Predicted values DB = Data used Pow0 Exp0 Exp1 BMA Predictive log score 34.109 35.244 33.968 31.394 Smallest The smaller the less information is lost

15. 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

16. GEOSTATISTICAL CHARACTERIZATION OF HIERARCHICAL MEDIA

17. Hierarchical Media Sedimentary deposits Fractured rocks Ritzi et al. 2006 Barton 1995

18. Hierarchical Media Properties are Multiscale (Nonstationary) Fractal plots for fracture trace maps (Barton 1995)

19. Permeability of Upper Oceanic Crust(Fisher 2005)

20. Longitudinal Apparent (Fickian) Dispersivity(Neuman 1995)

21. Hierarchical Media Properties appear to fit TraditionalStationary 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.

22. 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 = stationaryvariogram of nonstationary fractal field sampled on given support scale over finite window

23. Fractal Decomposition Power variogram is weighted integral of exponential or Gaussian modes: mode number integral scale

24. Fractal Decomposition Introducing upper and lower cutoffs: Lower cutoff renders field statistically homogeneous. Sample variogramsfittingstandard models of homogeneous fields often fit equally well.

25. TPVs are hard to distinguish from traditional stationarymodels Exponential (solid) compared with G-TPV (dashed) G-TPV = TPV based on Gaussian modes H = Hurst scaling coefficient

26. Tübingen Case Example 5 pumping tests in B wells interpreted graphically (Neuman et al. 2004, 2007): Integral Scale(ln T) = 2.5 m Variance(ln T) = 1.5 Support = 0Window = 25 m Note IS/W = 1/10!

27. Tübingen Case Example 312 flowmeterK values from all B and F wells were locally upscaled to obtain Variance(ln T) = 1.4 Window = 306 m Insufficientdata to estimate integral scale

28. 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

29. 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!): Y estimates Estimation variance

30. 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: Y estimates Estimation variance

31. Anomalous Transport in Randomly Heterogeneous Media

32. Contaminant Plumes in Groundwater at Hanford (2005)

33. Classic Fickian Modelfor Passive Tracer • Dispersiveflux = = constants • Advection-dispersion equation (ADE)

34. Classic Fickian Model Plume is Gaussian Its variance grows linearly with time

35. 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

36. Fieldmanifestation ofnon-Fickiandispersion AtBordendispersivity of solute slug under mean-uniform flow varies: • Linear in early travel-time/mean-distance (1:1 log-log slope) • Longitudinalstabilizes at “Fickian” asymptote • Transversepeaks, 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

37. Laboratorymanifestation ofnon-Fickiandispersion 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:

38. Simulated(localized stnADE) Effects of Boundaries & Conditioning on D (Morales-Casique et al. 2005)

39. Field+labmanifestations ofnon-Fickiandispersion Preasymptotic scaling also observed in radial field test(Peaudecerf and Sauty 1978) and2Dlab tests (Silliman and Simpson 1987; Sternberg et al. 1996) Apparent dispersivitiesworldwideincrease with scale at supralinear rate (Neuman 1990, 1995) excluding calibrated models (diamonds)

40. Laboratorymanifestations ofnon-Fickiandispersion 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)

41. Pore-scalesimulation ofnon-Fickiandispersion(Zhang & Lv 2007)

42. Space-time nonlocal theorybased 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 • Hence ADE d =  - scale dispersion tensor g = random source is stochastic.

43. 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) ensemblemean of c = optimum unbiased predictor of c = deterministic = smooth relative to c = (conditional) ensemblevariance of c = measure of predictive uncertainty

44. stnADE Transport Equation V = (conditional) meanv Q = (conditional) mean dispersive flux C, V, G=averagesover MCrealizations C, V, G, Qdefined on scale

45. stnADE Transport Equation(2) • Q = integro-differential(nonlocal)  = (conditional) nonlocalparameter • stnADE providessecond moments (predictive covariance) ofCand solute flux

46. Computational Example

47. Transport Conditions • Boundary conditions • Initial concentration

48. MC (left) and moment equation (right) velocity statistics

49. 2nd order (left) and iterative 2+ order (right) C

50. Longitudinal Mean Mass Flux