Geostatistical Inversion Methods

1. Geostatistical Inversion Methods Peter K. Kitanidis Stanford University With many thanks to students, post-docs, and other collaborators.

2. What this talk is about: • Motivation for down-scaling • Inverse problems • Geostatistical / Bayesian approaches • Quantification of uncertainty issues • Computational issues • Emphasis on core concepts

3. Premise Generally speaking, if you can afford to resolve processes at a finer scale, you have a better model. Resolve what? In space: Parameters, like conductivity, variables like pollutant concentration. In time: E.g., contaminant source distribution.

4. Site Characterization at a remediation location Courtesy of Murray Einarson How can one get the information to achieve resolution like that?

5. Measurements in Wells • Geologic Logging: At the surface, visual inspection and measurements of what comes out of the borehole during drilling and development. • Hydraulic and Geophysical logging: By instruments lowered into the borehole. Now, a plethora of methods are available. Challenge: Limited coverage!

6. Tomographic Methods • Measure integrated responses between “sources” and “receivers”. Feed data into computers. • Promise: Less intrusive than drilling many holes and yet gives spatially extensive information. • Challenges: • Limited resolution! • Often measuring not what we are interested in.

7. Practically all geophysical methods have been applied in hydrogeologic and environmental problems Courtesy of N. Linde

8. Insufficient Information • Estimate functions, like conductivity, from sparse and noisy observations; • The unknowns are sensitive to data gaps or flaws (Problem is ill-posed in the sense of Hadamard); • Usually, it is the smaller-scale variability that is hard to resolve.

9. Stochastic (aka Statistical or Probabilistic) Methods • Recognize that there are multiple solutions consistent with data. • Assign weights (“probabilities”) to these solutions.

10. Real (but unknown) Estimate A (fits data) Estimate B (plausible) Estimate C (a bit of both) Cheney, M. (1997), Inverse boundary-value problems, American Scientist, 85: 448-455.

11. The StructureKey idea: In any specific application,we know a number of things about the structure of the unknown (e.g., hydraulic conductivity) … • Values should be within a range. • There should be some degree of “continuity” and a reasonable “topology”. • Etc.

12. Bayesian Inference Applied to Inverse Modeling Information from observations Combined information (data and structure) Information about structure y : measurements s : “unknown”

13. How do we get the structure pdf? • Use geostatistical tools. • Or, use the maximum-entropy formulation to translate statements about structure into pdfs; • Use an “empirical Bayes” or “hierarchical Bayes” in which the structure pdf is parameterized and inferred from the data. • E.g.: We use cross-validation to find the parameters of the variogram used.

14. 20/30 Front View of the Sandbox 20/30 Data courtesy of Walter Illman Liu, X., and P. K. Kitanidis (2011), Large-scale inverse modeling with an application in hydraulic tomography, Water Resour. Res., 47(2), W0250110.1029/2010wr009144.

15. Back View of the Sandbox

16. Data Fitness

17. Validation with an extra pumping test

18. Tomogram with Half Million Pixels

19. Inverse Modeling: There is a Revolution going on! • Sensors • Fewer or no mechanicals parts, more electronics, more fiber optics, more telemetering. • Computers • They started becoming “useful” for inversions and tomography in the mid eighties. • Computing • High performance computing is the third pillar.

20. High Performance Computing • Exploiting the architecture of modern computers (supercomputers and computer clusters) to solve computationally large problems; • Exploiting structure in mathematical objects. E.g., a tridiagonal matrix is much easier to compute with. • So, adapt formulations and use the right numerical methods. • Most critical feature: How does computational effort scale with number of observations n and number of voxelsm, at high numbers? • n3 and m2 (even worse m3) are lethal; strive for n and m.

21. Faster Algorithm Strategies • Formulations that produce sparse matrices; • Iterative methods and approximate solutions; • “Covariance” and “Fisher Information” structure exploitation; • Ensemble strategies; • Data compression (like temporal moments or anchor points); • For nonlinear problems, more efficient iterative methods (reducing the number of forward runs).

22. Fast Computational Linear Algebra(with Darve, Ambikasaran, and Saibaba) • We are adapting to the inverse problem new numerical methods with near linear scaling using fast linear algebra. • Next, we will review the core ideas regarding the fast multipole method (FMM).

23. Core Ideas in FMM • Matrix-vectormultiplication, , is expensive, orderm2, when m is very large and A is full. • If A were low rank, then or visually • thenthe matrix vector multiplication can be done in order m operations. • Key idea: Exploit low-rank approximation of sub-blocks of the matrix. This reduces computational cost and potentially storage cost.

24. Cross-well, Straight-Ray Tomography with 10 Sources and 10 Receivers and a Very Fine Discretization Ambikasaran’sWork

25. Geostatistical Inversion For each source-receiver pair we have an observation yi, i= 1 to 100. The unknowns s are 10000x10000 = 108 y= Hs with Q the prior covariance (108 x 108) of s, X (108 x 1) the drift matrix of s, and R (102 x 102) the covariance of observation error. The system to be solved is:

26. Synthetic “True but Unknown” Estimated from data Direct and “exact” implementation of equations, exploiting only sparsity of H. The computation of HQHT took 183 seconds. Estimated from data Implementation of FMM tailored to this problem. The computation of HQHT took 6 seconds.

27. Hydraulic Tomography Work of Cardiff and Saibaba

28. Challenge: • The relation between observations (head) and unknowns (hydraulic conductivity) is nonlinear. The cost of computing the Jacobian matrix can be prohibitive . • In transient tomography, the number of observations is very large, though each measurement may contribute very little. We are trying alternative solution strategies

29. Transient 2-D Hydraulic Tomography log(conductivity) True Estimated Preliminary results to test alternative solver: Single well in the center, 25 observation wells on the boundary, practically no noise. Saibaba’s work

30. Total Variation (TV) Prior in Geostatistical Approach • Prior information : a few abrupt changes ins • θ₁,θ₂ are estimated using EM approach Lee’s work

31. 1-D Source Identification Problem

32. Identification of Zone Structure • Linear Tomography • Zone structure with intra-zone variability

33. Other methods Level Set Method Linear Variogram

34. Performance Comparison • * Thresholding is used for Geostatistical and TV approaches

35. Concluding Remarks • Inverse Modeling, an old field, is more exciting than ever. • We now can deal with many more data and more sophisticated models. • We are working on developing better algorithms to solve very large stochastic problems.

36. Supplementary Material

37. For example … We consider that the expected value of a norm of s is given; Then the pdf that maximizes the Shannon entropy functional can be written as

38. Or, consider another norm; Then the pdf that maximizes the Shannon entropy functional can be written as

39. The Likelihood Function (1/2) The relation between y and s is represented through: where e is the deviation between observation and prediction and is usually (but not necessarily) modeled as Gaussian,

40. The Likelihood Function (2/2) The likelihood is

41. Structural Parameters • θ1 is a structural parameter (or hyper-parameter related to structure); • θ2, variance of the measurement error, is another structural parameter.

42. Assume tentatively that structural parameters are known and we use one quadratic norm. Apply Bayes theorem. Find the s that minimizes this expression to determine the maximum a posteriori probability (MAP) solution.

43. Quantification of Uncertainty • Approximate, linearized analysis to get the posterior covariance. • Approximately, the ith conditional realization can be obtained through minimization of where and s’i and e’iare realizations according to their probabilistic models (prior pdf and measurement equation).

44. Bayesian Inference Including Structural Parameters y : measurement, s: unknown, and θ: structural-parameter vectors One approach: Find the q that maximizes the a posteriori probability. Finding the q from y is an algebraically overdetermined problem.