A Bayesian Approach to Earthquake Source Studies

1. A Bayesian Approach to Earthquake Source Studies Sarah Minson Caltech LLNL February 22, 2011

2. The people who make this possible • Mark Simons (Caltech) • James Beck (Caltech)

3. The big picture questions (or what I won’t be presenting) • Things we want to know • Do regions of co-seismic slip overlap with areas of post-seismic or inter-seismic slip? • How do hypocenter locations relate to co-/post-/inter-seismic slip? • Do earthquakes have smooth slip distributions and short slip durations? Vice versa? • Do earthquakes rupture at super-shear velocities? The answers to these questions come from finite fault earthquake source models…

4. Delouis et al. (2009) 2007 Mw 7.7Tocopilla,Chile Seismic + Static Teleseismic Strong motion Joint km Loveless et al. (2010)

5. Motivation • Seismology: • To get to earthquake physics, we need better source models • Math: • Many geophysical problems are under-determined • Ill-posed inverse problems require regularization • Roughness of the slip distribution (or whatever quantity is being regularized) cannot be identified because it was set a piori. What if you didn’t have to evaluate the inverse problem?

6. Bayes’ Theorem (1763) • For inverse problems: Posterior Prior Data Likelihood

7. Advantages of Bayesian analysis

8. Problems • Calculating posterior PDF generally requires Monte Carlo simulation • “Curse of Dimensionality” • Huge numbers of samples required for high-dimensional problems • Sampling can be inefficient especially in high-dimensional problems

9. Solution • Efficient parallel sampling algorithm • Metropolis algorithm and Markov Chain Monte Carlo (MCMC) are serial • Must efficiently sample highly anti-correlated model parameters • Must share information between worker CPUs

10. Cascading Adaptive Tempered Metropolis In Parallel: CATMIP • Tempering (A.K.A. Simulated Annealing)* • Dynamic cooling schedule** • Resampling** • Parallel Metropolis • Simulation adapts to model covariance** • Simulation adapts to rejection rate*** • Cascading * Marinari and Parisi (1992) ** Ching and Chen (2007) *** Matt Muto

11. CATMIP • Sample P(θ) • Calculate β • Resample • Metropolis algorithm in parallel • Collect final samples • Go back to Step 2, lather, rinse, and repeat until cooling is achieved

12. CATMIP Target • Sample P(θ) • Calculate β • Resample • Metropolis algorithm in parallel • Collect final samples • Go back to Step 2, lather, rinse, and repeat until cooling is achieved

13. Target CATMIP

14. Marginal PDF based on 400,000 samples of a 10-dimensional mixture of Gaussians

15. The seismic model • For each patch on a regularly gridded fault plane • Two components of slip • Rise time • Rupture velocity • Also may solve for ramp components of InSAR • Total number of parameters: • 4*Npatches (+ ramp)

16. Prior distribution • Slip • Rotated coordinate system relative to teleseismic rake direction • Slip parallel to rake has a positivity constraint • Forbidden from back-slipping more than 1 m • Prior distribution on slip perpendicular to rake is a zero-mean Gaussian • Prior distribution is populated from Dirichlet distribution • Prior distribution on rupture velocity and rise time is uniform

17. Forward model for static data • dpredicted=G*slip

18. Forward model for kinematic data • Green’s functions are pre-computed • These are then convolved with all possible source time durations and stored in memory • Convolution is too slow to be done at evaluation time • Final predicted waveforms are linear combination of each point source (with the appropriate pre-convolved source-time functions) scaled by slip, time-shifted according to results of Fast Sweeping • Fast Sweeping Algorithm (Zhao, 2005) is used to calculate initial rupture time at each source from Vr on each patch • Source parameters from each patch are interpolated onto a fine grid of point sources

19. Computation • Run parameters: • Markov chains = 500,000 • Steps per Markov chain = 100 • 2001 cores on CITerra Beowulf cluster • Results: • Cooling steps = 39 • Forward model evaluations = 1.931 billion • Run time = 10.6 hours • Total computation time = 76.3 million CPU seconds

20. How parallel is it? • Total computation time = 76.3 million CPU seconds • Each model evaluation = 0.035 sec • Static forward model = 0.005 sec • Kinematic forward model = 0.03 sec • Total model evaluation time = 68 million CPU seconds • ~90% of time spent on parallel model evaluation

21. In the future… • Caltech Center for Advanced Computing Research (CACR) is working on a GPU version of CATMIP • Michael Aivazis • Martin Michaelson • As of now, they have a first cut version with partial utilization of GPU • ~35x speed-up relative to CPU version • Need to move more of code onto GPU

22. Mw 7.7 Tocopilla, Chile Earthquake

23. The Data • Static GPS displacements • 1 Hz GPS time series • 6 interferograms

24. Static Prior Static Posterior/ Kinematic Prior Kinematic Posterior Cascading

25. Static vs. Joint Slip Distributions

26. Rise Time vs. Rupture Velocity

27. InSAR fits

28. GPS fits

29. Source parameters (now with uncertainties)

30. SourceProperties

31. Peak Slip: 3.015 ± 0.2733 m

32. Slip Heterogeneity

33. Visualization • How do you represent an N-dimensional PDF? • People like looking at slip models… • …But individual models can be very misleading…

34. Confidence bounds

35. Individual models (bad)

36. Even the mean and median tend to accentuate roughness…

37. Summary • CATMIP algorithm allows the sampling of very high-dimensional problems • Also useful for low-dimension problems with expensive forward models • Wide variety of uses in geophysics • Fully Bayesian finite fault earthquake source modeling • Resolution of the slip distribution and rupture propagation • Uncertainties on derived source properties • Determine which source characteristics are constrained and which are not

38. Future work • Data errors and model prediction errors are important • Assumed data errors control posterior model errors • Posterior distribution is also affected by model prediction errors (the failings of the forward model) • The next step is to estimate these errors

39. Thank you!