Scientific Questions
The outcomes & infrastructure arising from this work are relevant to answering several priority science questions in the Earth Sciences (National Academies, 2020), especially on plate tectonics, earthquakes & geohazards.
I. What is the relative strength of mantle heterogeneity at various spatial scales, & what does it imply for dynamic flow & dominant processes globally & regionally beneath North America?
Hypotheses - Velocity variations have the strongest power at the longest spatial scales, both globally (Figure 2) & regionally beneath the US. Dynamic processes in the deep mantle have a strong large-scale influence on most Cenozoic tectonic & magmatic activity in the western US. Temperature anomalies inferred in the western US are systematically under-predicted by regional models that do not capture heterogeneity at the longest spatial scales.
Global VS Heterogeneity spectra. Note dominance (red) of degrees 1-6 in the uppermost mantle & of degree 2 in transition zone & lowermost mantle.
II. Where are tomographic ‘blind spots’ that need to be plugged with future seismic deployments? Have theoretical improvements circumvented data coverage limitations? Relatedly, do the emerging & legacy tomographic techniques that use processed versus derived measurements “see” the same heterogeneity?
Hypotheses - Poorly sampled oceanic basins & regions with weak heterogeneity (e.g. passive margins) are blindspots. Theoretical improvements provide greater refinement on heterogeneity amplitudes than on patterns, but limitations due to station coverage persist. Structural models from quality-controlled derived measurements (e.g. Moulik et al., 2022) can be used as starting models in inversions of processed data (e.g. Full Spectrum Tomography, see Moulik and Ekström, 2014, 2016). Large-scale variations down to a nominal wavelength remain unchanged with legacy or emerging methods & afford similar fits to data.
Processed waveform (Left) & spectrum (Right) at Charters Towers, Australia from MW 9.1 2004 Sumatra Earthquake. Derived data types are listed; those in yellow are array-based techniques. Do tomographic techniques that use processed versus derived measurements “see” the same heterogeneity?
III. Is there a way to integrate existing knowledge on global & regional models? How well do multiscale models fit all types of data?
Hypotheses - Large-scale variations from global models are needed as baseline structure in multiscale models. Short-wavelength patterns can be integrated from regional models, but their amplitudes need to be adjusted through data validation. Using a robust three-dimensional Earth model (REM3D ) as the reference model in regional studies makes heterogeneity amplitudes & related interpretations more accurate, facilitating easier incorporation into the evolving ATLAS3D suite.
Scatter density plot of raw propagation phase anomaly measurements for 100 s R1 waves common to two techniques. Scatter points are colored based on high (red) or low (blue) spatial density. Strong agreement is found; minor discrepancies are from discrepant theoretical approximations, reference Earth models & processing schemes (Source: Moulik et al., 2022).
Background and Progress
The main goal of this project is creating the scalable infrastructure that assimilates robust structural features across a range of spatial scales into an evolving multiscale model suite (ATLAS3D), which represents the evolving state of our knowledge on the physical properties in the Earth’s interior. This work builds on several recent advancements in global tomography led by Moulik with co-workers worldwide.
Full Spectrum Tomography & Reference Earth models
Progress in tomographic modeling is driven by diverse data, ranging from astronomic-geodetic constraints to full seismic waveforms & derivative measurements of body waves (~1 – 20s), surface waves (~20 – 300s) & normal modes (~250 – 3000s). Full-spectrum tomography (FST) employs these observations to constrain physical properties – seismic velocity, anisotropy, density, attenuation & topography of discontinuities - in variable spatial resolution (Moulik & Ekström, 2014, 2016). Joint inversions account for different depth sensitivities to physical parameters & also reduce inter-parameter tradeoffs. Moulik has pioneered the FST approach for over a decade, & its utility is evidenced by its use as the de facto standard for inverting next-generation reference Earth Models for the community. Baseline large-scale structure of ATLAS3D is REM3D, which is derived using FST & reference seismological datasets.
Full Spectrum Tomography (FST) models the full spectrum of seismic data spanning 1-3000 s to constrain different types of heterogeneity in variable spatial resolution. Average Earth data, marked by *, ensure that bulk properties of 3D models are optimized simultaneously.
Bayesian Inference & Machine Learning with Rapid Forward Solvers
Our approach to ATLAS3D involves physics-informed machine learning (PIML) & Bayesian inference (posterior∝ likelihood x prior), where all information is represented in probabilistic terms (Box & Tiao 1973; Smith 1991; Tarantola, 2005). Residuals (R) are calculated between predicted (d) & observed data (dobs), either through a linear (G) or non-linear operator (g), which encapsulates the theoretical relationships that relate model parameters (m) to measurable quantities (Box 1). We have formulated a sparse matrix implementation for solving linear forward operators (NM, SW, TT), benchmarked on Graphical Processing Units (GPUs) for rapid evaluations of model likelihoods. Gaussianity & independence is assumed for (subsets of) observations & L2‐norm is used as misfit (𝝌2), resulting in diagonal data noise covariance matrix CD. The a priori knowledge on model m is incorporated through correlation between parameters (VP ⇔VS) in the matrix CM (Moulik & Ekström, 2016), & using REM3D as prior model of baseline structure (mprior). Sampling the full posterior distribution through random walks is prohibitive for many model parameters, even with dimensionality reduction of model vector m though basis functions. We instead extract features from existing Earth models & calculate joint likelihoods for the adversarial step of model selection. Various misfit measures (variance reduction, AIC, BIC) are tracked to reject scenarios. A limitation of this guided approach is that model selection of ATLAS3D is inherently focused on maximizing likelihood, biased towards existing knowledge, & underestimates true uncertainty of model parameters.
Formulae used for evaluating ATLAS3D model scenarios using Bayesian inference
Reference datasets in collaboration with the community
Reconciliation of data for the community REM3D project involves retrieving missing metadata, archiving in scalable storage formats, documenting outliers indicative of the limitations in some techniques, & quantifying summary reference data with uncertainties. Moulik et al. (2022) constructed a reference data set of multimode dispersion measurements by reconciling large & diverse catalogs of Love-wave (49.65 million) & Rayleigh-wave dispersion (177.66 million) from eight research groups worldwide. Reference data set summarizes measurements of dispersion at periods between 25-250 s of fundamental mode surface waves & up to six overtone branches from 44871 earthquakes recorded on 12222 globally distributed stations. These reference data, & others for body waves & normal modes, are used for model validation.
Scatter density plot of raw propagation phase anomaly measurements for 100 s R1 waves common to two techniques. Scatter points are colored based on high (red) or low (blue) spatial density. Strong agreement is found; minor discrepancies are from discrepant theoretical approximations, reference Earth models & processing schemes (Source: Moulik et al., 2022).
Utilizing legacy teleseismic & array data
Complexity in the wavefield also arises from the refraction & scattering of waves from heterogeneity (e.g., Masters et al., 1984; Lay & Kanamori, 1985). Directionality of the wave or its arrival angle, measured as deviation from the great-circle path, is sensitive to the gradients in velocity structure (Woodhouse & Wong, 1986). Figure below shows arrival angle deviations on the USArray Transportable Array (TA) from two source regions using the mini-array method, expanded after Foster et al. (2014). Arrival angle anomalies can arise at the boundary or outside the array (e.g. coastline interactions), but heterogeneity inside the array can also modify values along the path (e.g. Yellowstone). Details of these banded patterns can in turn constrain the amplitude & patterns of heterogeneity. Owing to their stronger sensitivity to small-scale structure than phase measurements, arrival-angle deviations provide a different perspective on wave propagation towards validating the ATLAS3D suite.
Arrival-angle anomaly for Rayleigh waves at 50 s period, for events located near Tonga (left) & Mexico (right). Black lines are the direction of wave propagation along the great-circle path. Banded patterns can arise outside the array & be modified by heterogeneity within, providing constraints on small-scale structure.
New wavelet-based method
We formulate new methods for feature extraction & assimilation. Figure below shows a cubed sphere with cartesian wavelet transform, modified from Simons et al. (2009). Advancements include: (1) Matrix implementation of transformations between geographic & wavelet domains with basis functions, (2) Windowing to account for extent of features, modulated with thresholding, & (3) Characterizing regions sans large spectral leakage.
Features of our multiscale wavelet configuration (four-tap Daubechies D4, resolution parameter N=8, scale J=5) for describing lateral variations (a) Cubed sphere centered on North America showing edges of 6 faces in red with 6 X 22N elements, (b) Approximate extent of spatial scales (1.5,3,6,12 deg.) overlain on topography, (c) Normalized power distribution with spherical harmonic degree, (d) Patterns & normalized amplitudes of orthogonal sets of scaling function, representing baseline structure, & wavelets for smaller-scale variations. Note that dashed circles for scales 2 & 1 (1.5 & 3 deg.) represent smaller variations than those for other scales.
Limits of wave propagation approximations
Waveform validation requires efficient evaluation of physical properties on the Gauss–Lobatto–Legendre (GLL) points of the global spectral-element mesh. These can be populated by evaluating the basis functions (spherical & cubic B-splines) of a starting model, permitting one-to-one mapping of locations to basis coefficients. Recent iterations of the adjoint tomography (GLAD-M25; Lei et al., 2020) have reported smaller-scale heterogeneity than those in the starting model (S362ANI, Kustowski et al., 2009). Meanwhile, revision to S362ANI has progressed towards including normal-mode splitting measurements (S362ANI+M; Moulik & Ekström, 2014), joint inversions of other parameters such as VP & density (Moulik & Ekström, 2016), & including reference datasets towards REM3D construction (Moulik et al., 2022). Figure below is a simulation from an earthquake in Argentina recorded at the Black Forest Observatory in Germany, with surface-wave windows used in GLAD-M25 construction shaded in blue. We modified SPECFEM3D_GLOBE code to facilitate comparisons of waveform fits from S362ANI+M & GLAD-M25, showing visual similarity, likely due to the dominant influence of large-scale structure. ATLAS3D waveform validation helps ascertain limits where structural deviations from REM3D perceptibly affect the measured wave packets. This procedure leverages the processing knowledge contained in the windows used for adjoint tomography.
Observed seismograms (black) & those predicted from large-scale heterogeneity (degree < 18, S362ANI+M, red) & GLAD-M25, which contains additional finer-scale features (blue). Maximum displacement (microns), misfits (M), correlations (C) & scaling (S) are provided to the right of each component trace.