Representation of Global Earth Models

1. Representation of Global Earth Models Interoperability, efficiency and prediction accuracy Bill Menke Lamont-Doherty Earth Observatory

2. Interoperability Being able to use someone else’s model in a quantitative way

3. Goalproviding the right information to make your model interoperable

4. even simple, well known issuescan cause substantialinteroperability problemshere’s an example

5. the world is not a sphere …. to center down north pole, Z f q X (X,Y,Z) Earth-Centered, Earth Fixed Cartesian coordinates Geographical latitude f not equal to geocentric latitude q Down not the same as direction to center

6. Recommendation – ECEF Coordinates You should always specify know to translate from what ever coordinate system that you’re using to earth-centered, earth- fixed (ECEF) Cartesian coordinates, (X,Y,Z) (0,0,0) at center of mass [0,0,1] points geographic north [1,0,0] points to intersection of equator and prime meridian And you should describe how to do it for any models that you publish or make widely available

7. Example: formula for converting geodetic coordinates to EDEF

8. ellipticity correction ~0.5sover continental distances

9. issues associated with ellipticity • Broad consensus to use geodetic lat, lon’s But what datum (model for earth’s shape)? GPS satellite system: WGS-84 hundreds of others in use Might effect: station locations event locations registration of regional tomographic results with rest of world

10. positioning errors of a few hundred meters are possible from using wrong datum significant in earthquake location; somewhat significant in global tomography

11. issues associated with ellipticity 2. How is depth defined? A problem here … “geodetic height at a point is the distance from the reference ellipsoid to the point in a direction normal to the ellipsoid” so height is parallel to local up/down, which makes intuitive sense

12. Oops! … depth defined this way misses the center of the earth to center down center Yet if you define depth as straight-line distance to the earth’s center, then geographic (lat, lon) varies along that line.

13. issues associated with ellipticity 3. How so you “ellipsify” a radial model of earth structure? You need to know how ellipticity varies with depth. Kennett, B., Seismological Tables: AK135 “… ellipticity coefficients were made for the iasp91 model using the algorithms presented by Doornbos (1988) with the density distribution from PEMC model of Dziewonski et al (1975) and the assumption that the ellipticity is nearly hydrostatic” OK: So how do I figure out what ellipticity vs. depth function Kennett used? I need to know e(r)! Why different?

14. issues associated with ellipticity 4. What is meant by a traveltime observation? A. Traveltime B. Traveltime corrected for somebody’s model of ellipticity I suspect that the reference to iasp91 in the AK135 Tables means that traveltime data were de-ellipsified using before being used

15. recommendation Good Model prediction + correction = observed data Bad Model prediction = observed data – correction bad because it tempts you to publish “data” that have been “corrected” in some not-very-obvious way. corrections: ellipticity, station elevation, anelasticity …

16. Case Study #1 of Interoperabiltycomparing continental scale models of P and S velocityInteroperabilty includes having enough information to understand the comparison

17. P velocity: Phillips (2007) isotropic model database of traveltimes 0.5x0.5 grid, registered on 0, ½ 50 km depthS velocity: Nettles & Dziewonski (2007) Voigt average of anisotropic model Love & Rayleigh waveforms 0.5x0.5 grid, registered on ¼, ¾ 50 km depth Here’s the biggest problem …

18. My solution: interpolate S to 0, ½ registration using bicubic splines Irony is that I know Nettles & Dziewonski (2007) model was based on spherical splines, so exact values at desired registration were “theoretically” available. Error probably small, since grid appears to be oversampled, but its hard to be sure.

19. 1:3 slope

20. Clearly some systematicdisagreement in location western edge of cratonwhat would we need to know to understand it ?

21. there’s a ratherlot of scatter arounddVs/dVp3expected fromeffect of temperaturewhat would we need to know to understand it ?

22. Case Study #2 of Interoperabiltyhow well does surface wave derived shear velocity model predict S-wave traveltimes?

23. Shear velocity: Nettles & Dziewonski (2007) Voigt average of anisotropic model Love & Rayleigh waveforms Crust 2.0 crust included in inversion but not given 50 km depth intervals given form mantle 0.5x0.5 grid, registered on ¼, ¾Traveltime calculation 3D model with 0, ½  registration, 50 km depth slices raytracing by shooting in model with tetrahedral splines Crust2.0 crust on top of Nettles & Dziewonski (2007) mantleS-wave Traveltimes Database of North America Regional Earthquake Traveltimes provided by W-Y Kim (personal communication) and where did he get the corresponding event locations?

24. Map of ray exit points craton Western NA Arrivals from Transition zone

25. Another Oops … (but for a different reason) The surface wave model has negative lithospheric velocity gradient over most of North America So there are no S-wave rays turning in the lithosphere

26. This problem is not limited to Nettles and Dziewonski’s (2007) North America model. Kutowski’s (2007) Eurasian model has regions with strong negative velocity gradients, too.

27. The Nettles and Dziewonski (2007) model and body wave travel times are inherently interoperable Why, we don’t know … maybe Earth looks different at very long periods or S-arrivals are not real rays (Thybo & Perchuc, 1997) About the best that one can do is to use it to “motivate the construction” of a model that would be useful for body waves

28. Example of “motivating the construction” Nettles & Dziewonski (2007) Vs-voigt at 50 km depth, scaled to Vp, and used to tweak velocity up & down uniformly through-out an ak135 lithosphere Travetime residuals predicted by model Traveltime residuals observed for Gnome explosion

29. many different model representations spherical harmonics voxels splines, with varying arrangement of nodes interpolant etc.

30. How should models be convertedbetween representations? Guiding principle a model is a way of summarizing the data therefore the conversion should try to reproduce an “idealized” version of the data used to create the model “idealized in the sense that one might not know exactly what data were used to create the model

31. Example crust mantle The figure at the left shows a hypothetical “layer cake” crustal model. Many such models have been published on the basis of active-source seismic refraction surveys. Suppose that we want to switch to a representation that uses linear splines. We use as the guiding principle that the layer-cake model probably fits the first arrival traveltime data well, in the distance range of a typical refraction experiment (say 200 km).

32. crust mantle calculated from simple formula Traveltime Prediction Layer Cake Model

33. crust mantle grey – tt’s from linear spline model calculated from inversion solid – tt’s from layer-cake model Comparison of layer-cake and linear spline models Traveltime “observations” and predictions The fit to the “data” is excellent, r.m.s.=100 ms, even with the crust being represented with only two linear splines. The down-side is that an inverse problem needed to be solved to computer the new model representation.

34. recommendation Model representations always be converted to preserve “idealized data”, Even though this requires solving an inverse problem However, the in hard cases the sense of idealization can be broadened to include preserving quantities that are merely “data-like”, rather than exactly data.

35. Example – traveltime data traveltime-like data: integrals of slowness along a suite of ray paths, whose shape is itself determined by the model traveltime-like data: integrals of slowness along a suite of prescribed curves that look something like ray paths

36. efficiency and prediction accuracy Some of my experiences with a ray-based traveltime tomography - earthquake location algorithm that uses a 3D model representation based upon linear tetrahedral splines

37. Choice of linear splines and tetrahedra motivated by efficiency

38. Advantages of tetrahedral models Ray paths known function (arc of circle) within tetrahedron Important ray integrals can be performed analytically, e.g. traveltime, T Where v is velocity, g is its gradient and

39. Disadvantages of linear tetrahedral models Curved surface of approximated as surface of polyhedron (but ½  node spacing gives reasonably accurate – 100 ms - traveltime). Velocity gradient discontinuous between tetrahedra (makes geometrical amplitudes rather rough) Not obvious how to generalize method to anisotropic earth models

40. A few examples

41. Curved interface approximated as surface of polyhedron Low Velocity Zone (LVZ) Slice through a 3D model represented with linear tetrahedral splines

42. Some sample ray paths

43. Fan array in next slide rays loop in LVZ Moho triplication upper crustal triplication Map view of ray exit points for source at

44. Traveltime plot for fan array geometry traveltime, s LVZ reverberations delayed by LVZ Fan array Y-distance

45. Efficiency issues Given an irregular tetrahedral grid .. 1. How do you figure out whether a given (x,y,z) point is in a given tetrahedron ? 2. How do you figure out which tetrahedron a given (x,y,z) point is in ?

46. 1. Is point in a given tetrahedron? Yes or No Its is if for each of the 4 faces of the tetrahedron it is on the same side of the face as the excluded vertex* of the face This one’s not One of the 4 faces Its excluded vertex * Three of the four vertices of the tetrahedron define a face. The fourth vertex is the “excluded” vertex. This point is on the same side as the excluded vertex

47. The inside/outside test needs to be very efficient, because it is used so often Thus information needed to perform it (e.g. the outward pointing normal of each face) needs to be pre-computed and stored.

48. 2. Which tetrahedron is a given (x,y,z) point in ? The probability is overwhelmingly high that its in the same tetrahedron as the last point that you focused upon … … because most operations are spatially localized An if not, then its probably in a tetrahedron that is close-by

49. Finding the tetrahedron by walking toward it current point test point New point Always move to tetrahedron adjacent to face with outward pointing normal is most parallel to the line connecting test point and new point vector connecting old and new points

50. Deciding which tetrahedron is adjacent to a given face of a given tetrahedron needs to be very efficient, because it is used so often Thus adjacency (nearest 4 neighbors) information needed needs to be pre-computed and stored.