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

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

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

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

### 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 (

**d**), either through a linear (

_{obs}**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 &

*L*‐norm is used as misfit (𝝌

_{2}^{2}), resulting in diagonal data noise covariance matrix

**The**

*C*_{D. }*a priori*knowledge on model

**m**is incorporated through correlation between parameters (V

_{P }⇔V

_{S}) in the matrix

**(Moulik & Ekström, 2016), & using REM3D as prior model of baseline structure (**

*C*_{M}**m**). Sampling the full posterior distribution through random walks is prohibitive for many model parameters, even with

_{prior}*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.

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

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

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

### 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 V_{P} & 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.