GlobalTomo: A global dataset for physics-ML seismic wavefield modeling and FWI

## Abstract

Abstract Global seismic tomography, taking advantage of seismic waves from natural earthquakes, provides essential insights into the earth’s internal dynamics. Advanced \ac fwi techniques, whose aim is to meticulously interpret every detail in seismograms, confront formidable computational demands in forward modeling and adjoint simulations on a global scale. Recent advancements in \ac ml offer a transformative potential for accelerating the computational efficiency of \ac fwi and extending its applicability to larger scales. This work presents the first 3D global synthetic dataset tailored for seismic wavefield modeling and full-waveform tomography, referred to as the \ac dataset dataset. This dataset is uniquely comprehensive, incorporating explicit wave physics and robust geophysical parameterization at realistic global scales, generated through state-of-the-art forward simulations optimized for 3D global wavefield calculations. Through extensive analysis and the establishment of \ac ml baselines, we illustrate that \ac ml approaches are particularly suitable for global \ac fwi, overcoming its limitations with rapid forward modeling and flexible inversion strategies. This work represents a cross-disciplinary effort to enhance our understanding of the earth’s interior through physics- \ac ml modeling.

## 1 Introduction

Global seismic tomography is a crucial yet intricate field within the earth sciences that facilitates understanding the earth’s internal structure and dynamics. This discipline has wide-ranging applications, from the discovery of natural resources (Virieux and Operto,, ) and the assessment of seismic hazards (Zang et al.,, ) to the exploration of our planet’s evolutionary history (Dziewonski and Romanowicz,, ). It utilizes seismic signals–earthquake-induced ground vibrations–recorded at surface stations to invert and map the earth’s interior. The fundamental challenge in this field lies in the integration of forward modeling of seismic wave propagation with \acfwi, both of which are essential for accurately correlating seismic data with velocity structures (Romanowicz,, ; Fichtner,, ). However, accurately simulating seismic waves through the Earth’s complex structures and inverting this data is challenging due to the heterogeneous nature of terrestrial media and the inherently ill-posed characteristics of \acfwi (Tromp,, ; Lyu et al.,, ).

Modern numerical simulation methods, such as finite difference methods (Kelly et al.,, ), finite element methods (Zienkiewicz et al.,, ; Moczo et al.,, ), and spectral element methods (; ; Krischer et al.,, ), are fundamental in the forward modeling of seismic waves. Although these methods have proven effective, they impose significant computational demands, particularly when deploying higher-resolution models necessary for detailed geophysical interpretation. High-resolution imaging requires the simulation of high-frequency waves, which in turn necessitates smaller time steps and finer spatial discretization. This is essential to prevent numerical dispersion and maintain the stability of the simulations, as mandated by the Courant-Friedrichs-Lewy condition (De Moura and Kubrusly,, ). Consequently, as the desired resolution increases, so does the computational burden, exponentially escalating the required computational resources and time, thereby introducing a substantial bottleneck ().

The computational intensity of forward modeling significantly constrains the efficiency of \acfwi. \acfwi aims to iteratively refine a model of the earth’s structure by minimizing the discrepancies between observed and simulated seismic data (Plessix,, ). The prolonged computation times required for forward modeling restrict the exploration of the model parameter space and limit the number of feasible iterations (Bozdağ et al.,, ; Lei et al.,, ). This limitation is particularly problematic when attempting high-resolution inversions over extensive areas. Consequently, this situation highlights the critical need for more efficient methodologies in both seismic forward modeling (Nissen-Meyer et al.,, ) and inversion strategies (Adourian et al.,, ).

Recent advances in \acml have shown significant promise in mitigating the aforementioned computational challenges in seismic tomography. Techniques such as \acpinn (Cuomo et al.,, ) and Neural Operator Learning (Lu et al.,, ; ; Wang et al.,, ) efficiently model complex physical systems by approximating solutions to the underlying \acppde. These methods have been successfully applied in fields including fluid dynamics prediction (Cai et al.,, ), new material discovery (Vasudevan et al.,, ), climate change modeling (O’Gorman and Dwyer,, ; Pathak et al.,, ), and enhancing biomedical treatments (Liu et al.,, ). By significantly reducing computational time compared to traditional numerical simulations, these \acml techniques enable more efficient high-resolution seismic modeling and inversion. Furthermore, modern \acml approaches can generalize beyond trained scenarios, facilitating the exploration of previously uncharted Earth’s interior.

However, the efficacy of \acml methods in geophysical applications is heavily contingent upon the availability of extensive, high-quality training data. Currently, there is a significant lack of benchmark synthetic datasets specifically designed for global seismic tomography and the \acml community. Although datasets to support \acfwi (Deng et al.,, ; Feng et al.,, ) have been introduced, they are mostly limited to subsurface exploration scenarios and do not fully represent the Earth’s entire scale, which is crucial for high-resolution wave propagation modeling. This deficiency highlights a critical gap in the resources needed to advance \acml applications in this area.

In response to these challenges, we introduce \acdataset, the first comprehensive 3D global full-waveform dataset specifically crafted for the \acml community. \acdataset encompasses high-resolution seismic simulations that span from the earth’s surface to its core. It features three distinct data tiers, integrating both acoustic and elastic wave equations to cater to a broad range of seismic phenomena. The dataset’s scope varies from a local scale of a 1-km radius simulated at 20Hz to a global scale encompassing the entire earth’s radius of 6,371 km at a 30-second period. We employ spherical harmonics to effectively parameterize velocity structures and account for variable seismic sources. Extensive benchmark analyses demonstrate \acdataset’s capabilities in fostering generalization, accelerating forward modeling, and addressing the inherent ill-posedness of seismic inversion tasks. By providing \acdataset, we aim to bridge the existing gap between geophysical sciences and the \acml community, thereby enabling more precise and efficient seismic modeling and inversion processes through \acml-driven innovations (Mousavi and Beroza,, ).

### 1.1 Related Work

##### \acl fwi

fwi has been extensively applied across a variety of domains, including shallow crustal imaging (Tape et al.,, ; Górszczyk et al.,, ), continental tomography (Chen et al.,, ; Zhu et al.,, ), medical imaging (Schreiman et al.,, ; Wiskin et al.,, ), and nondestructive testing (Rao et al.,, ). In seismology, classical \acfwi employs numerical simulations to derive seismograms from an initial earth model and iteratively optimizes this model using the adjoint-state method (Tarantola,, ; Tromp et al.,, ). While it has achieved significant success in both local and global inversions (Fichtner et al.,, ; French and Romanowicz,, ), \acfwi still confronts substantial challenges such as high computational demands and susceptibility to local minima (Tromp,, ). Emerging \acml approaches, bolstered by comprehensive global seismic datasets, present promising avenues to address these challenges, potentially leading to more efficient and precise \acfwi solutions.

##### Neural Operator Learning

Neural operator learning has garnered significant attention for solving parametric \acppde, as evidenced by a range of pioneering studies (Kovachki et al.,, ; ; ; Wang et al.,, ; Lu et al.,, ). Among these, DeepONet stands out for its flexibility in modeling mesh-free spaces and producing continuous solutions (Lu et al.,, ; Wang et al.,, ). Conversely, the \acfno operates with mesh-based inputs and outputs yet efficiently models the frequency domain and demonstrates rapid convergence during training (; Lu et al.,, ). Once trained, neural operators can swiftly predict solutions for new parameter sets with just a single forward pass. This capability is particularly advantageous in seismic applications, enabling the instant modeling of wave propagation under various velocity structures during inference. Such efficiency significantly streamlines the inversion process when compared to \acpinn methods (Song and Alkhalifah,, ; Rasht-Behesht et al.,, ), offering a more direct and rapid approach to seismic data interpretation.

## 2 Dataset Construction

### 2.1 Problem Definition

The primary goal of forward modeling here is to solve the 3D acoustic and elastic wave equations governed by interior \acppde and boundary conditions detailed in \crefsec:constraints. The objective is to predict the wavefield $\phi_{c,s}(p,t)$ at any spatial location $p=(x,y,z)$ and time $t$ across various velocity structures $c$ and source configurations $s$. This requires constructing and learning a nonlinear operator $G$, which maps velocity and source information into a wavefield output function, $G(c,s)$. The velocity structure $c$ is captured by sensor measurements at predetermined points, expressed as $c=[c(p_{1}),c(p_{2}),...,c(p_{m})]$. The source $s$ includes spatial locations and moment tensors. Given inputs $p$ and $t$, the function $G(c,s)$ computes the wavefield $\phi_{c,s}(p,t)$. This approach models wave propagation dynamics, aiding in accurate simulations and deeper geophysical understanding.

In particular, we monitor the displacements received by surface stations, known as seismograms, to infer the earth’s underlying structures. The inversion problem can be formulated as:

$$ $c^{*}=\argmin_{c}\|\mathbf{d}-\Phi(c,s)\|^{2}+\lambda F(c,s).$ (1) $$

Here, $\mathbf{d}$ represents the observed seismogram data. The function $\Phi(c,s)$ computes the seismograms based on the velocity structure $c$ and the source $s$. The term $F(c,s)$ is a regularization function that integrates prior knowledge and assumptions about the structure, while $\lambda$ is a coefficient that modulates the influence of the regularization term. The notation $\|\cdot\|^{2}$ denotes the L2 norm, used here to measure the misfit between the observed data and the predictions.

Traditional \acfwi relies on numerical simulations to compute $\Phi(c,s)$ and uses adjoint methods for optimizing $c$. Due to the intensive computational demand of these simulations, the optimization process is typically limited to a small number of iterations. We view \acml methods as a potent tool to expedite the forward simulation process and enable more flexible optimization strategies. The solutions of modern \acml approaches in addressing inversion problems is further discussed in \crefsec:pinn_nol.

### 2.2 Model Configuration

For generating global wavefield and seismogram data for our dataset, we employ AxiSEM3D (Leng et al.,, ), a forward simulator renowned for its efficiency in simulating global seismic wave propagation. AxiSEM3D is exceptionally adaptable, allowing customization to various model complexities (Leng et al.,, ). This simulator effectively utilizes the axisymmetry of global wavefields in a source-centered frame and extends calculations into the azimuthal Fourier domain. Further details are discussed in \crefsec:azimuth_expansion.

Our dataset comprises three tiers of increasing complexity, each supporting different training scales and applications. The first tier simulates acoustic wave propagation through a 1-km radius fluid sphere using a 20Hz mesh, providing a foundation in basic wave physics. The second tier extends these simulations to elastic wave propagation through an isotropic solid medium of the same scale, with source variations. This scale is comparable to 600 km continental simulations at a 30-second period, allowing exploration of complex P/S converted interactions. The third tier targets real-world earth applications, integrating acoustic and elastic simulations from the earth’s surface to its core over a 30-second period, addressing planetary-scale wave propagation for detailed geophysical analyses.

Each model’s 3D structure within these tiers is developed by introducing perturbations ranging from -10% to 10% to a 1D background model. These perturbations are meticulously designed to mirror the subtle yet significant 3D heterogeneity of the earth’s interior, as uncovered by global tomographic studies (Dziewonski and Woodhouse,, ; French and Romanowicz,, ; Thrastarson et al.,, ). Such adjustments are crucial for demonstrating the impact of tomographic features on earth’s geochemical and geodynamical processes, including mantle convection (van der Hilst et al.,, ), subduction of tectonic plates (Sun et al.,, ), and deep mantle heterogeneity such as the large low-shear-velocity provinces (Cottaar and Lekic,, ). Further details on the configurations of each model are provided in \crefsec:config.

Perturbations within each defined layer are parameterized using real spherical harmonics (Dahlen and Tromp,, ) up to degree 8, with radial values interpolated linearly between layers. Although some current mantle tomographic models offer higher resolution (Ritsema et al.,, ), spectrum analysis indicates that the significant power predominantly resides at lower degrees (Meschede and Romanowicz,, ; Ritsema and Lekić,, ). Consequently, our choice to limit structural parameterization to degree 8 in the published dataset effectively captures the predominant long-wavelength heterogeneity of the earth’s interior (Su and Dziewonski,, ). This approach reduces the number of parameters and enhances the solvability of this inverse problem.

### 2.3 Data Generation

tab:dataset details simulation inputs and outputs. Using \aclhs, we generated a set of random velocity structures and source configurations for extensive parameter coverage. Simulations were performed on a Dell PowerEdge R740 cluster with Intel Xeon Scalable processors (2,336 cores, up to 2.9GHz, 13,056GB RAM). Forward modeling for Acoustic and Elastic tiers required about 2,800 and 8,400 CPU hours respectively, while seismic simulations for the Real Earth tier needed approximately 100,000 CPU hours. Despite the high computational demands, this method is more cost-effective than traditional adjoint inversions (Thrastarson et al.,, ).

**Table 1: Simulation input and output size across three dataset tiers. The upper section of the table provides basic information about the three data tiers, including specific parameters and configurations. The lower section of the table details the dimensions of time, spatial elements, and predicted seismic data for each tier, respectively.**
| \topruleCategory | Variables | Acoustic | Elastic | Real Earth |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| \midrule\multirow5*Input | Sample number | 10,000 | 30,000 | 10,000 |  |  |  |  |  |  |
|  | Background model | Uniform | Uniform | PREM |  |  |  |  |  |  |
|  | Time range | 0-3 s | 0-3 s | 0-6000 s |  |  |  |  |  |  |
|  | Source | 0 | 6 | 9 |  |  |  |  |  |  |
|  | Structure | 405 | 1215 | 5427 |  |  |  |  |  |  |
| \midruleDimension |  | Time | Element | Seis | Time | Element | Seis | Time | Element | Seis |
| \midrule\multirow4*Output | Wavefield | 15 | 58,368 | 1 | 15 | 58,368 | 3 | 20 | 263,520 | 3 |
|  | Fourier wavefield | 15 | 3,648 | 16 | 15 | 3,648 | 48 | 20 | 16,470 | 48 |
|  | Seismogram | 150 | 1,369 | 1 | 150 | 1,369 | 3 | 6,000 | 703 | 3 |
|  | Storage | 88.54 GB | 657.28 GB | 1795.24 GB |  |  |  |  |  |  |
| \bottomrule |  |  |  |  |  |  |  |  |  |  |

Our dataset comprises two types of output data. The first type is surface seismogram data, which includes one-dimensional time series recordings of ground displacement at each surface station. These seismograms are observable globally and are crucial for the inversion of the earth’s internal structure. The second type is seismic wavefield data, capturing the propagating seismic wavefield within the earth’s interior. Although this spatially dense wavefield cannot be directly measured and used for traditional inversion, it reveals complex interactions between seismic waves and the earth’s three-dimensional structure. The intricacies of the wavefield, which contain rich details, have recently gained attention in global wavefield studies (Leng et al.,, ) but are often too challenging to identify manually. We provide this intermediate wavefield data to facilitate the training of neural networks. For detailed configurations of each type of output data, please refer to \crefsec:config.

## 3 Experiments

We experiment with several \acml baselines as the general solution operator for the 3D acoustic and elastic wave equations to accompany \acdataset. These baselines, trained on the Acoustic and Elastic tiers, serve as practical workflows for forward modeling tasks, predicting wavefields and seismograms based on given velocity structures and source parameters. Wavefield predictions include 7 timesteps from 0 to 3 seconds, with each timestep consisting of 16 slices and 3,648 points per slice. Seismogram predictions cover a series of time series of length 150 for 1,369 stations. These forward modeling baselines are not only utilized to simulate data but also employed in inversion tasks to refine the velocity structures based on their predictions.

### 3.1 Baselines

We implement our baselines using NVIDIA Modulus Sym (Hennigh et al.,, ), a framework designed for high-performance scientific computing.

##### \ac mm

The \acmm baseline calculates the average results from the training data and uses these averages as predictions for the test data. Although this method does not capture complex data patterns, it serves as a useful benchmark for assessing the performance of more sophisticated models. Additionally, it provides insights into the overall distribution of the data.

##### \ac mlp

The \acmlp baseline consists of a 6-layer fully connected neural network, each layer having a hidden size of 500 and utilizing SiLU activation functions. This model takes velocity structures as inputs and predicts both wavefields and seismograms.

##### \ac hfn

Considering the tendency of neural networks to favor low-frequency solutions (Rahaman et al.,, ), we explore the use of a \achfn. This model combines the gating mechanisms of Highway Networks (Srivastava et al.,, ) with the frequency domain learning capabilities of Fourier Neural Networks (Tancik et al.,, ). The \achfn is designed to handle deep learning tasks requiring both complex layer-to-layer transformations and frequency-specific learning, making it potentially more adept at managing intricate patterns in seismic data.

##### DeepONet

Instead of directly mapping velocity structures to fixed output grids, this approach involves building a framework that models both wave propagation and seismogram generation on a continuous spatial-temporal domain. We utilize the DeepONet (Lu et al.,, ), where the trunk net processes inputs of positions and times, and the branch net handles velocity structures. The outputs from both nets are combined through a dot product before passing through a linear output layer to produce the final prediction. To improve training efficiency and reduce memory, we developed an optimized method for parallel multiplication of the trunk and branch net outputs; see \crefpara:eff_deeponet.

##### \ac pido

While DeepONet is trained using traditional supervised learning, which might limit its generalization in high-resolution domains, we explore additional capabilities with \acpido (Wang et al.,, ). This model is designed to enforce physically consistent representations by incorporating interior and free-surface constraints during training. These constraints are integrated by calculating the necessary derivatives within the \acppde using automatic differentiation.

### 3.2 Results

#### 3.2.1 Forward Modeling

##### Generalization to unseen structures

To investigate the ability of the baseline models to generalize across unseen seismic structures, we trained them using approximately 90% of the velocity structures available in our dataset. The remaining 10% of the structures were used as a test set to evaluate the models’ performance. Quantitative assessments were carried out using \acrl2 and \acr metrics.

**Table 2: Quantitative results of baselines. Performance of baseline models in forward modeling is evaluated on the test set using \acsrl2 and \acsr. Results include both mean and standard deviation values. The average inference time for each model is also presented.**
| \toprule\multirow2*Wave | \multirow2*Model | \multirow2*Representation | Wavefield | Seismogram | \multirow2*\makecellTime |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  |  |  |  |  |  |  | (ms) |
|  |  |  | \acsrl2 $\downarrow$ | \acsr $\uparrow$ | \acsrl2 $\downarrow$ | \acsr $\uparrow$ |  |
| \midrule\multirow4*Acoustic | MM | Vector | 0.495$\pm$0.048 | 0.871$\pm$0.026 | 0.597$\pm$0.069 | 0.802$\pm$0.053 | - |
|  | \acmlp | Vector | 0.356$\pm$0.035 | 0.937$\pm$0.012 | 0.446$\pm$0.051 | 0.896$\pm$0.026 | 1.771 |
|  | H-Fourier | Vector | 0.397$\pm$0.054 | 0.917$\pm$0.066 | 0.594$\pm$ 0.073 | 0.820$\pm$0.068 | 2.604 |
|  | DeepONet | Point | 0.368$\pm$0.050 | 0.927$\pm$0.066 | 0.503$\pm$0.061 | 0.862$\pm$0.067 | 2.898 |
| \midrule\multirow4*Elastic | MM | Vector | 1.000$\pm$0.001 | -0.006$\pm$0.113 | 1.000$\pm$0.001 | 0.001$\pm$0.116 | - |
|  | \acmlp | Vector | 0.617$\pm$0.038 | 0.790$\pm$0.027 | 0.534$\pm$0.068 | 0.848$\pm$0.039 | 1.773 |
|  | H-Fourier | Vector | 0.735$\pm$ 0.061 | 0.689$\pm$0.048 | 0.613$\pm$0.067 | 0.790$\pm$0.049 | 2.397 |
|  | DeepONet | Point | 0.682$\pm$0.035 | 0.734$\pm$0.031 | 0.565$\pm$0.065 | 0.826$\pm$0.043 | 2.775 |
| \bottomrule |  |  |  |  |  |  |  |

Figure: Figure 2: Qualitative forward modeling prediction by \acmlp in the Acoustic tier. A single slice is displayed. The source is located at $x=0.00$ and $z=0.80$. The background illustrates the velocity structure. The seismogram depicts the time series received by stations around the surface of this slice.
Refer to caption: x2.png

Figure: (a)
Refer to caption: x3.png

Differences in baseline performance highlight the strengths of various modeling approaches. As detailed in \creftab:for_res, \acmlp achieves high correlation coefficients and low \acrl2 values, demonstrating its efficacy in learning from vector-based representations. Observations from training indicate that the \achfn model reaches low training errors quickly, suggesting its proficiency in capturing frequency-based signals; however, it tends to overfit, resulting in limited generalization capabilities compared to the \acmlp. While vector-based methods excel with structured data, they face challenges in predicting arbitrary coordinates and time steps. In contrast, the point-based DeepONet model offers greater flexibility in generalizing across different coordinates and time steps, supporting mesh-free predictions. Nonetheless, it requires more computational resources during training to achieve comparable performance to the \acmlp. Further training details are available in \crefsec:models. Qualitative evaluations, illustrated in \creffig:for, corroborate these findings. The \acml baselines closely match the actual data across several time steps, capturing the complex dynamics of wave propagation on unseen structures. As time advances, emulation errors tend to increase, reflecting the growing complexity of wave dynamics over extended periods, as depicted in \creffig:for_a.

##### Complexity analysis

The performance of \acmm can serve as an indicator of the dataset’s complexity and distribution characteristics. Smaller error margins suggest more consistent patterns across velocity structures. The Elastic tier, evidenced by the poorer performance of \acmm, exhibits greater complexity than the Acoustic tier due to the intrinsic properties of the elastic wave equation and variations in source dynamics. The trend towards increased complexity is also corroborated by the performance metrics observed across different \acml baselines.

##### Forward modeling speed

A single round of forward simulation using a numerical solver requires 120 seconds with 24 CPU cores. In contrast, \acml approaches achieve forward modeling speeds typically ranging from 1 to 3 milliseconds on a single GPU, making them about 60,000 times faster than traditional methods. Note that the acceleration rate may vary with the degree of parallelism. While most computational costs for \acml seismic modeling are incurred during the generation of training data, once trained, the neural operator can predict all wavefields in a single inference step. This capability underscores the potential of \acml approaches for real-time seismic modeling and extensive forward modeling applications. For a comparative visualization, see \creffig:for_b.

##### Generalization to higher temporal resolution

To evaluate the flexibility of \acfwi in real-world scenarios through continuous wavefield prediction across the entire temporal domain, we tested the generalizability of a trained neural operator to higher resolution snapshots. Specifically, we utilized DeepONet, which accommodates mesh-free inputs. DeepONet was trained on a uniform velocity structure using acoustic waves at timesteps 1, 3, 5, and 7, and subsequently evaluated on a denser sequence from timesteps 1 to 8, within the same temporal range. The findings, illustrated in \creffig:for_c, show that DeepONet, being purely data-driven, initially struggled with adapting to this higher-resolution domain. To enhance its generalization capabilities, we incorporated physical constraints (detailed in \crefsec:constraints) into its training regimen. This integration of physical principles significantly improved DeepONet’s predictions on previously unseen timesteps, reducing its dependency on extensive labeled training data sets for developing a super-resolution model.

#### 3.2.2 Inversion

In this section, we explore how the generalization and acceleration capabilities of \acml methods in forward modeling can significantly improve the inversion process. We present three strategies to address the challenges inherent in traditional \acfwi and enhance the efficiency of inversion operations.

##### Gradient-based optimization

r0.42
The inversion optimization process. Performance improves through optimization across 200 iterations. Higher degrees capture increasingly shorter-wavelength structures, enhancing model fidelity.

Following traditional \acfwi methods, we initially set fixed parameters for established baselines and iteratively refine the velocity model to minimize the L1 misfit between the predicted and actual seismograms. Our optimization strategy employs 200 iterations of gradient descent using the LBFGS algorithm, with a learning rate of 0.08 and a history size of 30. A regularization term is integrated to promote stable convergence. We explore the effectiveness of gradient information from \acml forward models in aiding the inversion optimization process. For each test structure, five initial points are randomly selected, and optimization is performed using the best-performing forward model (\acmlp in the Acoustic setting). The inversion performance is assessed using \acr across various structural scales with 80, 245, and 405 free parameters corresponding to degrees 4, 6, and 8, respectively. As shown in \creffig:inv_iteration_square, the correlation between the inverted model and the ground truth strengthens with each iteration step, confirming that \acml-derived gradients can effectively guide the inversion process within 200 steps.

##### Sampling with multiple starting points

r0.4
Inversion with increasing starting points. Correlation between the inverted models and the ground truth strengthens as the number of initial starting points increases, indicating improved robustness and accuracy of the inversion process.

The initial conditions significantly influence the performance of inversion. To enhance inversion, we experimented with starting the optimization from multiple initial points. Specifically, we used \acml forward modeling to generate 1,000 random starting points for each structure in the test set using \aclhs. Each starting point underwent optimization with the chosen algorithm, and the best result was selected for analysis. As illustrated in \creffig:inv_start_square, performance consistently improves with an increase in the number of starting points across various structural scales. This suggests that using more starting models could further boost inversion efficacy.

In summary, both increasing the number of optimization steps and using more starting models significantly enhance inversion performance. This underscores the ability of \acml emulation to effectively tackle the challenges of ill-posedness and local minima—issues that are often problematic in traditional numerical solvers due to the high computational demands of extensive forward modeling. The use of \acml approaches allows for more comprehensive exploration of solution spaces, which is critical for achieving accurate and robust inversion results.

##### Direct inversion mapping

An alternative inversion strategy is to directly map observed seismograms to velocity structures. This method bypasses the conventional requirement for a forward \acml model, potentially addressing the optimization challenges inherent in traditional \acfwi. We evaluated this approach by training an \acmlp with 6 layers and 1000 hidden units, designed to input acoustic wave seismogram data and output velocity structures. Post-training, this model predicted reasonable velocity structures across various degrees in an unseen test set, achieving an average \acr of 0.826. Visualization of these results is provided in \creffig:inv. This direct inversion mapping approach offers promising insights into solving \acml-based inversion problems more efficiently.

Figure: Figure 4: Qualitative inversion results. The results compare the inverted structure of the earth’s interior, derived from acoustic data, with the ground truth across different spherical harmonic degrees. The color scale, transitioning from blue to red, indicates negative to positive perturbations respectively. For additional visualizations, please refer to \crefsec:vis.
Refer to caption: x8.png

#### 3.2.3 Discussion and Future Work

- [leftmargin=*]
- •
Uniqueness of \acdataset.
The universality and representativeness of the dataset are critical for ensuring robust performance when applying the trained network to real-world data. In the realm of earth science, we benefit from concentrating on a single, unique object–the Earth–from its surface to its core. Unlike datasets from localized explorational or geological sites, which can vary significantly, our dataset is crafted to encapsulate the global attributes of the earth’s interior. This focus enhances the dataset’s applicability as all simulations and models are framed within a consistent setup. Future work could expand the scope by straightforwardly incorporating more degrees of spherical harmonics to explore the Earth’s interior with greater detail.
- •
Broader interest of \acdataset.
Beyond application in earth, our dataset is also crucial for advancing planetary sciences, especially for celestial bodies like Mars and the Moon. Seismic data for these planets are rare and costly to obtain. A pre-trained model that simulates seismic activity on these bodies can significantly improve our understanding of their interior.
- •
Potential of \acml.
\acml presents a swift and powerful approach to addressing long-standing scientific challenges. Utilizing a finite set of simulated data that incorporates physics-inspired efficient representations, \acml methods can perform global forward modeling and inversion across various structures and scales with minimal time investment. The accuracy and efficiency of these methods could be further enhanced by increasing the dataset size and adopting more advanced models.

## 4 Conclusion

We hope \acdataset, a comprehensive 3D global full waveform dataset for physics-\acml research, would advance seismic forward modeling and inversion. Our analysis showcases the ability of modern \acml methods’ capabilities in seismic wavefield modeling, delineating their strengths and limitations. The potential for discoveries in understanding the earth’s internal processes grows as we bridge the gap between theoretical \acml innovations and practical geophysical applications.

##### Societal Impacts

A deeper understanding of the earth’s internal structure brings numerous societal benefits, including improved disaster preparedness and mitigation, safer urban planning and infrastructure development, more effective and sustainable management of natural resources, and environmental preservation. This research also fosters global scientific and technological collaboration, further amplifying its positive societal impacts.

1.1

## Acknowledgments and Disclosure of Funding

The authors would like to thank Miss Chen Zhen (BIGAI) for making the nice figures, Prof. Yizhou Wang (Peking University) for supporting the project, and NVIDIA for their valuable technical support for accelerated computing and artificial intelligence. S. Li, Z. Mu, S. Xin, and Y. Zhu are supported in part by the NSFC (62376009) and the Beijing Nova Program. Z. Li and X. Song are supported in part by the National Key R&D Program of China (2022YFF0800601), China Postdoctoral Science Foundation (2023M740043), and the Emerging Engineering Interdisciplinary–Young Scholars Project, Fundamental Research Funds for the Central Universities. K. Leng is supported in part by the EPSRC grant, Blueprinting for AI for Science at Exascale (BASE-II, EP/X019918/1).

## Contents of Supplementary Material

linkcolor=black
\startcontents[sections]
\printcontents[sections]l1
\hypersetuplinkcolor=red

## Appendix A Data and Code Access

We have released our dataset on our website \hrefhttps://global-tomo.github.io/https://global-tomo.github.io/. The code for data loading, training, and evaluation has been documented and published on \hrefhttps://github.com/lishiqianhugh/GlobalTomo/https://github.com/lishiqianhugh/GlobalTomo/. We will continuously update and maintain both the code and data to ensure they reflect the latest advancements.

## Appendix B Physical Representations

### B.1 Physical Constraints

We describe the physical constraints of the seismic wave propagation that can be used in the training process of neural networks. For the acoustic dataset tier, the wave propagation follows the acoustic wave equation in the fluid domain $\Omega_{f}$:

$$ $\frac{\partial^{2}_{t}\chi}{\kappa}=\nabla\cdot(\frac{\nabla\chi}{\rho}),\quad \textit{in}\quad\Omega_{f}.$ (A1) $$

The output displacement $\mathbf{u}=\nabla\chi/\rho$, $\chi$ denotes the wave potential, $\rho$ represents the fluid density, and $\kappa$ represents the bulk modulus. The acoustic dataset has a pressure-free boundary condition that can be stated by:

$$ $\chi=0,\quad\textit{on}\quad\partial\Omega_{f}.$ (A2) $$

The elastic wave propagation in the solid domain $\Omega_{s}$ follows:

$$ $\rho\partial^{2}_{t}\mathbf{u}-\nabla\cdot(\mathbf{C}:\nabla\mathbf{u})=0, \quad\textit{in}\quad\Omega_{s}.$ (A3) $$

where $\mathbf{C}$ represents the elasticity tensor. The elastic dataset has a traction-free surface boundary condition:

$$ $\hat{\mathbf{n}}\cdot(\mathbf{C}:\nabla\mathbf{u})=0,\quad\textit{on}\quad \partial\Omega_{s}.$ (A4) $$

The Real Earth tier is a combination of these two scenarios. The seismic wave propagation in the fluid outer core follows \crefeq:eom_fluid while in the resting solid places follows \crefeq:eom_solid. This tier still has a traction-free boundary condition same as \crefeq:BC_solid. The solid-fluid coupling is controlled by the displacement- and stress-continuity condition on the fluid-solid interface $\Sigma$:

$$ $\hat{\mathbf{n}}\cdot\mathbf{u}=\chi_{r}/\rho\quad,\quad\hat{\mathbf{n}}\cdot( \mathbf{C}:\nabla\mathbf{u})=\partial^{2}_{t}\chi\hat{\mathbf{n}}\quad\textit{ on}\quad\Sigma.$ (A5) $$

### B.2 Expansion to Azimuthal Fourier Domain

Describing the 3D earth is memory-consuming. Inspired by the smoothness and symmetry found in the physical properties of the earth, representing the wavefield in the Fourier domain can achieve a more efficient description. If the time-dependent equations of motion can be generalized in computational domain $\Omega$ as:

$$ $\mathscr{L}\mathbf{u}=\mathbf{f}$ (A6) $$

The displacement $\mathbf{u}$, operator $\mathscr{L}$, and external forces $\mathbf{f}$ will be characterized to the required azimuth resolution defined by $n_{u}$:

$$ $\mathbf{u}(r,\theta,\phi;t)=\sum_{|\alpha|\leq n_{u}}\mathbf{u}^{\alpha}(r, \theta;t)\exp(\text{i}\alpha\phi)$ (A7) $$

$$ $\mathscr{L}(r,\theta,\phi;t)=\sum_{|\beta|\leq n_{u}}\mathscr{L}^{\beta}(r, \theta;t)\exp(\text{i}\beta\phi)$ (A8) $$

$$ $\mathbf{f}(r,\theta,\phi;t)=\sum_{|\gamma|\leq n_{u}}\mathbf{f}^{\gamma}(r, \theta;t)\exp(\text{i}\gamma\phi)$ (A9) $$

Thus this method reduces the dimension of azimuth that is needed to solve a true 3D problem and saves huge computation costs while also ensuring the accuracy of the solution. If we apply the structural parameterization of spherical harmonics and determine the highest degree of the model, AxiSEM3D currently has the best theoretical efficiency for calculating their global wavefields, which only need 2$n_{u}$ times the computation of a 2D problem. In \acdataset, we provide the Fourier series of the wavefield for building more robust forward models with efficient representation in the future.

### B.3 Weak Formulation and Dimension-Reduced Form

The weak formulation of the elastic equation can be obtained by taking dot product with an arbitrary test function $\mathbf{w}$ and integrating by parts over the volume $\Omega_{s}$:

$$ $\int_{\Omega_{s}}\Big{(}\rho\partial^{2}_{t}\mathbf{u}\cdot\mathbf{w}+\nabla \mathbf{u}:\mathbf{C}:\nabla\mathbf{w}\Big{)}\,\textit{d}\mathbf{x}^{3}=\int_{ \Omega_{s}}\mathbf{f}\cdot\mathbf{w}\textit{d}\mathbf{x}^{3}$ (A10) $$

The stress-free boundary condition is thus naturally merged into the weak formulation of the constraint.

Its corresponding dimension-reduced form in a 2D domain $D$ can be written down as:

$$ $\split\sum_{|\alpha|\leq n_{u}}[\langle\rho^{-(\alpha+\beta)}\partial^{2}_{t} \mathbf{u}^{\alpha},\mathbf{w}^{\beta}\rangle_{D}+a_{D}^{-(\alpha+\beta)}( \mathbf{u}^{\alpha},\mathbf{w}^{\beta})]=\langle\mathbf{f}^{-\beta},\mathbf{w} ^{\beta}\rangle_{D}\\ \forall\mathbf{w}^{\beta}\,\textit{in}\,D,\quad\forall\beta\in\{-n_{u},...,n_{ u}\}$ (A11) $$

$$ $\langle\mathbf{u}^{\alpha},\mathbf{w}^{\beta}\rangle_{D}:=\int_{D}u_{i}^{ \alpha}w_{i}^{\beta}s\mathrm{d}s\mathrm{d}z$ (A12) $$

$$ $a_{D}^{\gamma}(\mathbf{u}^{\alpha},\mathbf{w}^{\beta}):=\int_{D}u^{\alpha}_{i, j}C^{\gamma}_{ijkl}w^{\beta}_{k,l}s\mathrm{d}s\mathrm{d}z$ (A13) $$

where the azimuthal expansion of $\mathbf{u}$, $\mathbf{f}$, $\rho$, and $C$ follows:

$$ $\mathbf{u}=\sum_{|\alpha|\leq n_{u}}\mathbf{u}^{\alpha}(s,z;t)\Psi^{\alpha}(\phi)$ (A14) $$

$$ $\mathbf{f}=\sum_{|\alpha|\leq n_{u}}\mathbf{f}^{\alpha}(s,z;t)\Psi^{\alpha}(\phi)$ (A15) $$

$$ $\rho=\sum_{|\alpha|\leq n_{u}}\rho^{\alpha}(s,z;t)\Psi^{\alpha}(\phi)$ (A16) $$

$$ $C_{ijkl}=\sum_{|\alpha|\leq n_{u}}C_{ijkl}^{\alpha}(s,z;t)\Psi^{\alpha}(\phi)$ (A17) $$

Here, $\Psi^{\alpha}(\phi)$, $\alpha\in\{-n_{u},...,n_{u}\}$ denote the azimuthal Fourier modes,

$$ $\Psi^{\alpha}(\phi):=\exp(\mathrm{i}\alpha\phi),\mathrm{i}=\sqrt{-1}.$ (A18) $$

## Appendix C Earth Configuration Details

### C.1 Concept Illustration

- [leftmargin=*]
- •
Source: Described by location and moment tensor, a source of seismic energy can be either natural or artificial.
- •
Velocity Structure: Refers to the distribution of seismic wave velocities within the earth’s interior, determined by the physical and chemical properties of the materials.
- •
Wavefield: The spatial pattern of wave propagation within the earth.
- •
Seismogram: The recordings of ground motions received by stations on the surface.
- •
Forward Modeling: The process of computing numerical solutions of the wave propagation from a defined model of the earth’s interior.
- •
Full Waveform Inversion (FWI): An advanced seismic inversion technique that leverages the entire seismic waveform signals to develop high-resolution models of the earth’s interior.

### C.2 Three Dataset Tiers

We introduce the details of the three tiers generated by AxiSEM3D. The simulator is licensed under the MIT license at \hrefhttps://github.com/AxiSEMunity/AxiSEM3Dhttps://github.com/AxiSEMunity/AxiSEM3D.

Our dataset includes three tiers of synthetic seismic data. The first tier models acoustic wave propagation through a fluid sphere with a 1-km radius using a 20 Hz mesh. The second tier models elastic wave propagation through an isotropic solid medium of the same scale, with variations in the source. The third tier is designed for real-world earth applications, integrating both acoustic and elastic simulations over a 30-second period.

We export the 3D velocity perturbations in Cartesian coordinates to facilitate rapid interpolation and retain the original spherical harmonics coefficients for use in network training. To ensure the adherence to physical constraints, we disable attenuation.

##### Acoustic Ball

The acoustic ball model is designed on a 1000 m radius ball with pure fluid medium and a mesh accurate to 20 Hz. This set of models has a pressure-free boundary condition. The wave propagation in the model obeys the acoustic wave equation. The medium only has the P-wave velocity (vp) attribute. The 1D background model has a uniform vp value of 1.0 km/s across the ball. The 3D superimposed velocity perturbation is parameterized by spherical harmonics up to degree 8, leading to 81 spherical parameters per layer, multiplied by 5 layers, resulting in 405 free parameters in total.

The source is located at a fixed 200 m depth at the north pole with a fixed fluid pressure source.

##### Elastic Ball

The elastic ball model is also designed on a 1000 m radius ball with pure elastic medium and a mesh accurate to 20 Hz. This set of models has a stress-free boundary condition. The wave propagation in the model obeys the elastic wave equation. The isotropic medium only has three independent attributes: P-wave velocity (vp), S-wave velocity (vs), and density (rho). The 1D background model has a uniform vp value of 1.5 km/s, vs value of 1.0 km/s, and rho of 1.0g/cm^3 across the ball. The 3D superimposed velocity perturbation is parameterized by spherical harmonics up to degree 8, leading to 81 spherical parameters per layer, multiplied by 5 layers and 3 medium attributes, resulting in 1215 free parameters in total.

The source is located at a fixed 200 m depth at the north pole with a variable moment tensor source mechanism. This introduces 6 extra free parameters in source terms.

##### Real Earth

This set of models is designed for the ultimate real-earth applications. We utilize the isotropic PREM (Dziewonski and Anderson,, ) background model and a global mesh with a period accurate to 30 seconds that encompasses various earth’s layers including the crust, mantle, outer core, and inner core, each characterized by their unique seismic attributes. The Real Earth tier has a stress-free boundary condition. The 3D velocity perturbations are imposed on the earth’s internal 23 discontinuities, resulting in a total of 5427 free parameters up to degree 8 across the whole earth.

The source is located at stochastic depths between 0 - 800 km at arbitrary geographical coordinates with a variable moment tensor source mechanism. This introduces 9 extra free parameters in source terms.

##### Surface Seismogram Output

The seismogram data contains a one-dimensional time series recording ground displacement at each surface station. Each trace is low-pass filtered to 30 seconds to avoid numerical noise. We place such stations evenly along latitude and longitude to represent a spatially abundant coverage. The actual station deployments limit the real seismic data coverage, however, this issue can be considered by dumping part of the data during the application. The seismogram data can be observed globally and be used in the inversion of the earth’s internal structure.

##### Seismic Wavefield Output

The seismic wavefields are saved in discrete 15 or 20 time snapshots, (total time = 3 s or 6,000 s), represented by 16 azimuthal slices with 3,648 or 16,470 (from surface to the core) element points in each. Besides, we also saved the azimuthal Fourier coefficients of the global wavefields for a continuous and compact representation in three dimensions. The wavefields and the seismograms are saved as HDF5 files for fast loading.

## Appendix D Baseline Model Details

### D.1 Evaluation Metrics

The \acrl2 metric is particularly valuable in scenarios where it’s crucial to evaluate the performance of predictive models across datasets of varying scales or inherent variability. By normalizing the Mean Squared Error (MSE) with the variance of the actual values, this metric effectively renders the error measurement scale-independent. The formulation is as follows:

$$ $\text{Relative}L_{2}\text{Loss}=\sqrt{\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^ {2}}{\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}}}$ (A19) $$

where $y_{i}$ and $\hat{y}_{i}$ denote the actual and predicted values of the $i$-th observation, respectively. $\overline{y}$ is the mean of all actual values in the dataset, calculated as $\frac{1}{n}\sum_{i=1}^{n}y_{i}$.

We also use the Person Correlation Coefficient (R) to measure the similarity between the prediction and ground truth. The \acr is given by:

$$ $R=\frac{\sum(y_{i}-\overline{y})(\hat{y}_{i}-\overline{\hat{y}})}{\sqrt{\sum(y _{i}-\overline{y})^{2}\sum(\hat{y}_{i}-\overline{\hat{y}})^{2}}}$ (A20) $$

### D.2 Efficient DeepONet Design

The original DeepONet takes paired velocity structure and spatial-temporal points for the trunk network and branch network. The output of DeepONet is calculated by combining the outputs from the branch and trunk networks, typically through a multiplication operation. Specifically, if the output from the trunk network is $a=[a_{1},a_{2},\dots,a_{d}]$ and the output from the branch network is $b=[b_{1},b_{2},\dots,b_{d}]$, then the predicted value $\hat{\phi}(p,t)$ at output position $p$ and time $t$ is given by:

$$ $\hat{\phi}(p,t)=\sum_{i=1}^{d}a_{i}\cdot b_{i}(p,t)$ (A21) $$

Here, $a_{i}$ and $b_{i}(p,t)$ are the features obtained from the branch and trunk networks, respectively, and their product followed by summation forms the estimate of the output function.

During training, this network is optimized to minimize the error between the true output $\phi(p,t)$ and the predicted output $\hat{\phi}(p.t)$, typically using the mean squared error as the loss function:

$$ $\text{Loss}=\frac{1}{n}\sum_{j=1}^{n}\left(\phi(p_{j},t_{j})-\hat{\phi}(p_{j}, t_{j})\right)^{2}$ (A22) $$

where $n$ is the number of output points.

However, this is memory-consuming since the velocity structures are duplicated for each point. Indeed, one velocity structure can be paired with a group of spatial and temporal points $P,T$ simultaneously. Here we design a more efficient DeepONet that can predict the whole 3D outputs just using one forward of branch net. The calculation of the output function now is given by:

$$ $\hat{\Phi}(P,T)=\sum_{i=1}^{d}a_{i}\cdot B_{i}(P,T)$ (A23) $$

where $B_{i}(P,T)$ indicates a vector of the i-th latent feature from a batch of points in the same velocity structure.
$\hat{\Phi}(P,T)$ indicates a group of predicted output in the same structure.
Suppose the input dimensions of the velocity structure and the spatial-temporal point are $m$ and 4, respectively, the memory usage in the modified DeepONet can be $\frac{4+m}{4+m/n}$ times smaller than the original one. The training speed is about 3 times faster.

### D.3 Hyperparameters

For reproducibility, we present the training configuration of baselines in \creftab:wf_model and \creftab:seis_model. We employ the SiLu activation function and use the Adam optimizer for parameter tuning. The initial learning rate is set at 0.0003, and it decreases exponentially by a factor of 0.95 every 1000 steps. Both the \acmlp and \achfn models consist of six layers with 500 neurons each. The DeepONet model features six layers of 600 neurons each, and includes trunk and branch projection layers, each with 1000 neurons. The final output is computed via a dot product of the latent trunk and branch vectors.

The \acmlp and \achfn models are trained on a single A800 SXM4 80GB GPU. In contrast, DeepONet requires 8 A800 SXM4 80GB GPUs for training due to the substantial increase in data size associated with point-based prediction.

**Table A1: Wavefield model details.**
| \topruleWave type | Model | Layer num | Hidden size | Training steps | Input size | Output size |
| --- | --- | --- | --- | --- | --- | --- |
| \midrule\multirow3*Acoustic | \acmlp | 6 | 500 | 10000 | 405 | 408576 |
|  | \achfn | 6 | 500 | 3000 | 405 | 408576 |
|  | DeepONet | 6 | 600 | 100000 | 405+3 | 1 |
| \midrule\multirow3*Elastic | \acmlp | 6 | 500 | 30000 | 1221 | 408576 |
|  | \achfn | 6 | 500 | 7000 | 1221 | 408576 |
|  | DeepONet | 6 | 600 | 40000 | 1221+3 | 1 |
| \bottomrule |  |  |  |  |  |  |

**Table A2: Seismogram model details.**
| \topruleWave type | Model | Layer num | Hidden size | Training steps | Input size | Output size |
| --- | --- | --- | --- | --- | --- | --- |
| \midrule\multirow3*Acoustic | \acmlp | 6 | 500 | 10000 | 405 | 205350 |
|  | \achfn | 6 | 500 | 3000 | 405 | 205350 |
|  | DeepONet | 6 | 600 | 100000 | 405+3 | 1 |
| \midrule\multirow3*Elastic | \acmlp | 6 | 500 | 10000 | 1221 | 205350 |
|  | \achfn | 6 | 500 | 4000 | 1221 | 205350 |
|  | DeepONet | 6 | 600 | 100000 | 1221+3 | 1 |
| \bottomrule |  |  |  |  |  |  |

## Appendix E Discussion

### E.1 Comparison of Modern ML-based Inversion

Here we discuss the difference between modern \acml-based inversion approaches. The two most dominant schools are PINN and neural operator learning.

We first recall the formulation of the inversion problem.

$$ $c^{*}=\argmin_{c}\|\mathbf{d}-\Phi(c,s)\|^{2}+\lambda F(c,s)$ (A24) $$

In comparing the inversion strategies of PINN and operator learning, we observe distinct paths of optimization. The PINN-based approach seeks to optimize the model parameters $c$ in conjunction with the forward modeling function $\Phi$ (Rasht-Behesht et al.,, ). Thus, the $c$ and $\Phi$ are closely coupled and the inversion of different structures requires training from scratch. This method limits the model’s ability to generalize across different velocity structures.

In contrast, the operator learning approach addresses this limitation by initially learning a universal operator to model forward problems across all $c$ (Yang et al.,, ). Subsequently, it fine-tunes the model parameters $c$ for each specific case using the pre-trained $\Phi$. Once trained, the neural operator can achieve inversion on any structures in inference time.

Our study adopts the latter approach to achieve fast and flexible inversion. We study both the sampling-based and optimization-based methods when using the pre-trained neural operator to solve the inversion problem.

### E.2 Is Predicting Wavefield Necessary for Inversion?

As shown in \crefsec:inversion, direct inverse mapping demonstrates better inversion performance than the optimization-based method. We want to discuss whether predicting wavefield is still necessary from the aspects of complexity, interpretability, and flexibility.

##### Complexity

Direct \acml inversion methods offer an efficient alternative to traditional optimization-based approaches by using neural networks to approximate the inverse mapping from data to subsurface properties directly. This method significantly cuts down on computational time and resources, simplifying the traditional iterative inversion process into a single step. In contrast, splitting the \acfwi process into two stages might increase complexity regarding the implementation and integration of different computational frameworks or software tools. Additionally, more stages could lead to greater cumulative errors.

##### Interpretability

Predicting the wavefield explicitly as part of the inversion process can enhance the interpretability of the results. In \acfwi, the aim is to reconstruct the subsurface properties using wave propagation data. Although direct inverse mapping can produce a result of the inferred structure, it doesn’t explicitly model the intermediate process and tends to lack transparency in its operations. By decomposing the inversion into two distinct phases—predicting the forward wavefield and then optimizing the subsurface model—researchers can trace the effects of specific inputs or changes in model parameters on the outcomes. This division aligns closely with traditional FWI methods where each step of the process can be examined and understood separately. Furthermore, explicitly modeling the wavefield allows for a more detailed examination of the physical accuracy and stability of the inversion process, enabling researchers to validate and refine the underlying physical models being used.

##### Flexibility

Predicting the wavefield explicitly provides opportunities to apply a broader range of optimization strategies tailored to specific challenges or characteristics of the data. For example, if the predicted wavefield is treated as a separable component, researchers can implement advanced gradient methods such as adaptive step-sizing or preconditioning. These techniques can address specific dynamics of wave propagation in different geological settings, potentially enhancing the accuracy and efficiency of the inversion. Preconditioning, in particular, can help tackle issues like the ill-posed nature of the inversion problem or the presence of sharp local minima by altering the optimization landscape.

### E.3 Connection with Other Fields

Similar to the challenges of interpreting the human mind in the field of intuitive physics, the task of inverting Earth’s structure involves comparable approaches and discussions.

Both seismic inversion and intuitive physics deal with inverse problems—inferring causes from effects. Seismic inversion back-calculates to infer the Earth’s interior from surface measurements (Tromp,, ). Similarly, intuitive physics tries to deduce the cognitive rules or mental models that people use to predict physical outcomes from their observed behavior in experimental settings (Kubricht et al.,, ).

Regarding methodologies, both fields heavily utilize computational models. In seismic studies, numerical methods and statistical models predict how geological features influence seismic waveforms (Tarantola,, ; Tromp et al.,, ). Similarly, in intuitive physics, physics engines simulate how humans might predict physical interactions in their surroundings (Battaglia et al.,, ; Li et al.,, ). Machine learning, particularly neural networks, plays a significant role in both areas. In seismic inversion, neural networks can be trained to recognize patterns in seismic data that correlate with specific subsurface structures (Yang et al.,, ). In intuitive physics, neural networks model human prediction and reasoning processes about physical laws, learning to mimic human judgment ().

Uncertainty management is crucial in both domains. In seismic inversion, uncertainties arise from limited sensor coverage and noise in data. In intuitive physics, uncertainties stem from the variability in human perception and cognitive biases. Both fields develop methodologies to handle these uncertainties, often through probabilistic models (Tromp,, ; Battaglia et al.,, ).

## Appendix F Datasheet

### F.1 Motivation

- •
For what purpose was the dataset created? \acdataset was created to promote the physics-\acml research on solving seismic full waveform modeling and earth structure inversion.
- •
Who created the dataset (e.g., which team, research group) and on behalf of which entity (e.g., company, institution, organization)? The dataset was created by a group of geophysical scientists and \acml researchers.
- •
Who funded the creation of the dataset? The authors would like to thank NVIDIA for their valuable technical support for accelerated computing and artificial intelligence. S. Li, Z. Mu, S. Xin, and Y. Zhu are supported in part by the NSFC (62376009) and the Beijing Nova Program. Z. Li and X. Song are supported in part by the National Key R&D Program of China (2022YFF0800601), China Postdoctoral Science Foundation (2023M740043), and the Emerging Engineering Interdisciplinary–Young Scholars Project, Fundamental Research Funds for the Central Universities. K. Leng is supported in part by the EPSRC grant, Blueprinting for AI for Science at Exascale (BASE-II, EP/X019918/1).

### F.2 Composition

- •
What do the instances that comprise the dataset
represent (e.g., documents, photos, people, countries)? \acdataset contains many simulated wavefield sequences and seismogram time series documented as arrays in HDF5 files.
- •
How many instances are there in total (of each type, if appropriate)? The acoustic tier has 10,000 cases of simulation results. The elastic tier has 30,000 cases and the real elastic tier has 10,000 cases.
- •
Does the dataset contain all possible instances or is it
a sample (not necessarily random) of instances from a larger set?
Yes, the dataset contains all possible instances for training and evaluation.
- •
What data does each instance consist of? Each instance consists of inputs and outputs. The inputs have velocity structures and source parameters. The outputs contains three types of features including the wavefield, the Fourier wavefield, and the seismogram. The wavefield includes a sequence of snapshots and each snapshot consists of 16 slices, together with the 3D coordinates and timesteps. The Fourier wavefield resembles the wavefield but presents each 3D snapshot in terms of 16 Fourier series. The seismogram consists of time series data from multiple stations, which also provide the 3D coordinates of each station and corresponding timesteps.
- •
Is there a label or target associated with each
instance? Yes, each instance has clearly defined input and output variables.
- •
Is any information missing from individual instances?
No.
- •
Are there recommended data splits (e.g., training,
development/validation, testing)? We recommend using 90% of the data for training/validation and the remaining 10% for testing.
- •
Are there any errors, sources of noise, or redundancies
in the dataset? The velocity structures are replicated across various features to facilitate data loading and training. Given the small size of these structures, the redundancy in storage can be considered negligible.
- •
Is the dataset self-contained, or does it link to or
otherwise rely on external resources (e.g., websites, tweets,
other datasets)? The dataset is self-contained.
- •
Does the dataset contain data that might be considered
confidential No.
- •
Does the dataset contain data that, if viewed directly,
might be offensive, insulting, threatening, or might otherwise
cause anxiety? No.

### F.3 Collection Process

- •
How was the data associated with each instance
acquired? Was the data directly observable (e.g., raw text, movie
ratings), reported by subjects (e.g., survey responses), or
indirectly inferred/derived from other data (e.g., part-of-speech
tags, model-based guesses for age or language)? The data was collected from a set of simulations using a 3D global seismic model named AxiSEM3D.
- •
What mechanisms or procedures were used to collect the
data (e.g., hardware apparatuses or sensors, manual human
curation, software programs, software APIs)? We utilized a cluster of CPUs to run the AxiSEM3D.
- •
Who was involved in the data collection process (e.g.,
students, crowdworkers, contractors) and how were they compensated
(e.g., how much were crowdworkers paid)? Geophysics researchers in the author list were involved in the data collection process and no crowdworkers were involved.

### F.4 Preprocessing/cleaning/labeling

- •
Was any preprocessing/cleaning/labeling of the data done
(e.g., discretization or bucketing, tokenization, part-of-speech
tagging, SIFT feature extraction, removal of instances, processing
of missing values)? We examine the failure cases in simulation and remove them from training and evaluation. The data is processed into HDF5 files with labeled variables for easy loading.
- •
Was the “raw” data saved in addition to the preprocessed/cleaned/labeled data (e.g., to support unanticipated future uses)? Yes, the raw data is saved in our cluster.
- •
Is the software that was used to preprocess/clean/label the data available? Yes, we uploaded the code on GitHub.

### F.5 Uses

- •
Has the dataset been used for any tasks already? No.
- •
What (other) tasks could the dataset be used for? The dataset could be used for seismic forward modeling, earth structure inversion, and source inversion.
- •
Is there anything about the composition of the dataset or the way it was collected and preprocessed/cleaned/labeled that might impact future uses? No.
- •
Are there tasks for which the dataset should not be used? No.

### F.6 Distribution

- •
Will the dataset be distributed to third parties outside of the entity (e.g., company, institution, organization) on behalf of which the dataset was created? The dataset is open to the public.
- •
How will the dataset be distributed (e.g., tarball on website, API, GitHub)? The dataset will be distributed on Google Drive and the code will be published on GitHub.
- •
Will the dataset be distributed under a copyright or other intellectual property (IP) license, and/or under applicable terms of use (ToU)? The dataset will be distributed under the CC BY-NC-SA license.
- •
Have any third parties imposed IP-based or other restrictions on the data associated with the instances? No.
- •
Do any export controls or other regulatory restrictions apply to the dataset or to individual instances? No.

### F.7 Maintenance

- •
Who will be supporting/hosting/maintaining the dataset? The \acdataset group will support, host, and maintain the dataset.
- •
How can the owner/curator/manager of the dataset be contacted (e.g., email address)? The owner/curator/manager of the dataset can be contacted through the following emails: Shiqian Li (shiqianli@stu.pku.edu.cn) and Zhi Li (zl382@pku.edu.cn).
- •
Is there an erratum? No.
- •
Will the dataset be updated (e.g., to correct labeling
errors, add new instances, delete instances)? Yes, the dataset will be continuously updated if there is a necessity to improve accuracy and other related information. The updates will be released on the website.
- •
If the dataset relates to people, are there applicable
limits on the retention of the data associated with the instances
(e.g., were the individuals in question told that their data would
be retained for a fixed period of time and then deleted)? The dataset is not related to people.
- •
Will older versions of the dataset continue to be
supported/hosted/maintained? Yes, older versions of the dataset will continue to be supported, maintained, and hosted.
- •
If others want to extend/augment/build on/contribute to
the dataset, is there a mechanism for them to do so? Yes, the contributor can contact us through email.

## Appendix G More Visualization

### G.1 Forward Modeling

We further present the visualization of forward modeling on four novel velocity structures to demonstrate that the \acml baselines are capable of discerning differences in wave propagation. We display snapshots at timesteps 5, 7, 9, 11, and 13, illustrating that the variance in wavefields grows over time. Refer to \creffig:for_2,fig:for_3,fig:for_4,fig:for_5,fig:for_6 for detailed visualizations.

Figure: Figure A1: Snapshot 5 of forward modeling. We show the ground truth and predicted wavefield snapshot within four velocity structures: A, B, C, and D.
Refer to caption: x9.png

Figure: Figure A2: Snapshot 7 of forward modeling. We show the ground truth and predicted wavefield snapshot within four velocity structures: A, B, C, and D.
Refer to caption: x10.png

Figure: Figure A3: Snapshot 9 of forward modeling. We show the ground truth and predicted wavefield snapshot within four velocity structures: A, B, C, and D.
Refer to caption: x11.png

Figure: Figure A4: Snapshot 11 of forward modeling. We show the ground truth and predicted wavefield snapshot within four velocity structures: A, B, C, and D.
Refer to caption: x12.png

Figure: Figure A5: Snapshot 13 of forward modeling. We show the ground truth and predicted wavefield snapshot within four velocity structures: A, B, C, and D.
Refer to caption: x13.png

### G.2 Inversion

We present additional visualization of inversion results using various unseen seismogram data. The figures display both the actual and the inverted velocity perturbations across five depths, along with their correlations. Refer to \creffig:inv_198,fig:inv_350,fig:inv_578,fig:inv_600,fig:inv_638,fig:inv_662,fig:inv_744,fig:inv_831,fig:inv_865,fig:inv_909.

Figure: Figure A6: Example 1 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x14.png

Figure: Figure A7: Example 2 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x15.png

Figure: Figure A8: Example 3 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x16.png

Figure: Figure A9: Example 4 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x17.png

Figure: Figure A10: Example 5 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x18.png

Figure: Figure A11: Example 6 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x19.png

Figure: Figure A12: Example 7 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x20.png

Figure: Figure A13: Example 8 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x21.png

Figure: Figure A14: Example 9 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x22.png

Figure: Figure A15: Example 10 of acoustic inversion. The colors range from blue to red, representing velocity perturbations from negative to positive. There are five rows, each depicting the velocity structure at different depths: 0m, 200m, 400m, 600m, and 800m. The first column displays the actual velocity structure, the second column shows the structure as derived from direct inverse mapping, and the third column illustrates the correlation between the spherical harmonics of the actual structure and those of the inverted structure.
Refer to caption: x23.png

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. Could the authors clarify the intended benefit of including the full wavefield in the dataset, given that only surface measurements are typically available in realistic FWI scenarios? How does the inclusion of the full wavefield impact the dataset’s applicability to practical, real-world inversion problems where the wavefield is not accessible?


2. Could the authors provide more details on the distribution of samples within each tier of the dataset? Are there visualizations available to illustrate the diversity of scenarios represented?


3. Were any generalization tests conducted using direct inversion methods to assess the dataset’s performance across varied inversion scenarios? If not, it would be much appreciated if some generation tests could be conducted using the direction inversion methods used in this current manuscripts (InversionNet-3D, InversionMLP).
soundness: 2
strengths: This paper addresses a key gap by providing a global-scale dataset that integrates ML-friendly features with robust physical modeling for seismic applications. The dataset’s design across three tiers (acoustic, elastic, and real Earth) allows it to support scalable and complex seismic modeling tasks.
summary: The paper, "GlobalTomo: A Global Dataset for Physics-ML Seismic Wavefield Modeling and FWI," introduces a 3D global synthetic dataset designed specifically for machine learning applications in seismic wavefield modeling and full-waveform inversion. The GlobalTomo dataset integrates high-resolution seismic simulations that span various scales, from acoustic and elastic wave propagation to planetary-scale simulations representing real Earth conditions. The paper demonstrates the utility of this dataset by benchmarking several machine learning models for both forward modeling and inversion tasks, underscoring the potential of ML to accelerate seismic data interpretation and enhance our understanding of Earth’s interior.
weaknesses: 1. Wavefield Estimation. The necessity of including the wavefield as part of the dataset is not entirely convincing. While the authors provide further discussion on this in Supplementary Section F.2, this inclusion could be misleading for two primary reasons.

First, from a practical inverse problem perspective, only surface wavefield measurements (seismograms) are typically available in real-world FWI applications. The full wavefield throughout the Earth’s interior is generally inaccessible and, therefore, impractical for use in inversion.

Second, traditional model-based FWI methods do require calculating the full wavefield as they rely on directly solving the wave equation. However, end-to-end ML-based inversion strategies approximate the inversion process and thus do not depend on the full wavefield. These methods instead focus on estimating target parameters directly from surface measurements, bypassing the need for a complete wavefield solution.


2. Comprehensive Assessment of the Dataset.  An essential purpose of a benchmark dataset is to support various realistic scenarios that may arise in practical applications. A key example would be testing generalization capabilities. However, this aspect is not clearly addressed in the manuscript.

It remains unclear whether the dataset, or its three tiers, adequately represents different distributions or if they cover a range of scenarios that would facilitate generalization. No visualizations are provided to illustrate sample distributions, making it difficult to assess the dataset’s diversity. Additionally, I did not find any generalization tests conducted with direct inversion methods, as outlined in Section 3.2.2. This raises concerns about the dataset's robustness across diverse inversion tasks.

---

confidence: 5
contribution: 4
initial_rating: 6
presentation: 4
questions: - If you are sure all methods are well-trained and under the best hyperparameters, could you explain why MLPs perform better than methods like CNN and Transformer in 3D cases, which seems counterintuitive? Or, could you provide additional analysis or ablation studies that might shed light on why MLPs are outperforming other architectures?

- Could you compare the performance of these methods in scenarios of out-of-distribution?

- In line 396, the paper says “To evaluate the flexibility of FWI”. Why is this  FWI? I think this section is talking about forward modeling. Could you clarify whether this section is indeed about forward modeling rather than FWI?
soundness: 2
strengths: - The dataset itself is good and could benefit the community a lot. 
- The paper is well-written and easy to follow.
- The figures and tables are clear.
summary: - This paper presents the first 3D global synthetic dataset tailored for seismic wavefield modeling and full-waveform tomography, referred to as the Global Tomography (GlobalTomo) dataset.
- This paper also utilizes and compares several machine learning models on forward modeling and FWI based on this dataset.
weaknesses: I think the machine learning models are kind of simple. Based on my experience, MLPs seem not to be the best choices for 3D problems compared with other methods you chose in the paper (CNN or transformer). But the results show that MLP / InversionMLP performs best in most cases. I am concerned whether all methods are well-trained and under the best hyperparameters.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: 1.	Section 1 Introduction (Line 106): While the paper proposes to address the gap for high-resolution wave propagation modeling for the Earth’s entire scale, the Real Earth dataset employs wave with 30 second period. How does the 30 second wave period (likely low frequency) facilitate in generating high-resolution imaging of the internal Earth’s structure?

2.	Section 2.2 Model Configuration (Line 191): What is meant by Model Configuration (machine learning model, wave velocity distribution, or something else)? Assuming it means wave velocity distribution. The paper proposes to generate synthetic 3D velocity structure by perturbing a 1D background velocity model by -10 to 10%. While this may be more accurate for local heterogeneity, this assumption seems to oversimplify the 3D structure for the Earth’s scale data. Is there any rationale behind using this simplification or limitations in representing Earth’s complex structure?

3.	Section 2.3 Data Generation (Line 213): How many examples for velocity structures were generated per tier of dataset? How many source configurations were generated per tier of dataset? Is it same as number of velocities generated or different?

4.	Section 2.3 Data Generation (Line 221): The definition of the input and output variables is not clear and difficult to understand. For example, what does Source, Structure, Sample Number mean here? How does these parameters affect the quality of generated dataset in terms of heterogeneity, geological structures, and scale?

5.	Section 3 Experiments (Line 254): For wavefield dataset, there are 7 timesteps considered (0 to 3 seconds), with each timesteps consist of 16 slices and 3648 points per slice. How does these 16 slices represent the entire internal Earth’s structure? Does each slice represent a unique sub-region inside the Earth? (A schematic figure to represent this would be good) Why are there 3648 points per slice, are these points indicating positive and negative wave displacement in that slice? 

6.	Section 3.2.1 Forward Modeling (Line 324): Based on the wavefield dataset description in the line 253, the total number of timesteps should be 7 from 0 to 3 seconds, the Figure 2 shows that the results for 13 timesteps. Does the dataset contain 7 or 13 timesteps?

7.	Section 3.2.1 Forward Modeling (Line 324): Figure 2 shows seismogram data, how are source and receivers positioned relative to each other? A schematic diagram to explain this may be helpful to understand the seismogram data for each tier of dataset. 

8.	Section 3.2.1 Forward Modeling (Line 402): How was the physics constraints enforced to DeepONet that led to improved generalizability? Appendix B.1 shows governing physics equations for each tier of dataset but do not provide insights into training the model with physics. 

9.	Section 3.2.2 Inversion (Line 424) and Figure 4: What does it mean to select five random points for each test structure? Does this mean that there are 5 initial velocity guesses, and inversion is carried out for all of them independently? Which tier dataset is used for inversion here? Why the model pre-trained on acoustic tier is selected? Which model from the acoustic tier is selected for the forward simulation? Can we have some figures to compare the inverted velocity structures from different models for each dataset? How does the model pre-trained on one dataset generalizes to other dataset for both forward and inverse problems?

10.	Section 3.2.2 Inversion (Line 446): In the dataset description, the acoustic data is defined as a 1 km radius fluid sphere. Figure 6 shows inverted Earth’s structure derived from acoustic data. What does it mean? 

11.	Section 3.2.2 Inversion (Line 457): Does different starting point undergo optimization independently? If not, how the inversion algorithm is processing multiple starting point to reach to one ground truth while affecting each other?
soundness: 2
strengths: 1.	The idea to generate a large-scale dataset for global seismic tomography is interesting and essential to bridge the gap between machine learning and seismic tomography.
2.	The proposed dataset has three tiers to model seismic wave propagation from simple to complex settings and for each tier the dataset contains seismic wavefield and seismogram information.  
3.	The paper explores multiple ML methods to establish a benchmark for both forward and inverse problems in seismic tomography.
summary: The paper proposes datasets for the forward and inverse modeling for global tomography. The forward modeling entails predicting seismogram from a velocity structure, while the inverse modeling solves for a velocity structure from an observed seismogram. They have considered three different cases for data generation: (1) Acoustic Ball (1km radius sphere filled with pure fluid medium), (2) Elastic Ball (1km radius sphere isotropic media), and (3) Real Earth (based on isotropic PREM model). Each dataset has two types of output data, the first type is surface seismogram data, and the second type is the seismic wavefield data. In the paper, the authors have employed several ML models for forward and inverse problems and have shown the efficacy of ML models in solving the two problems on the proposed dataset.
weaknesses: 1.	Although the paper is overall well-written, there are certain sections which is little obscure like model configuration, spherical harmonics, and data generation. These sections need improvement to understand how the dataset was generated, model training and inference results.
(a)	For generating different examples, do we always perturb the same 1D background model or different? How is the data generation algorithm create different varying geological settings in the dataset? How do we measure the variability in different velocity structures?
(b)	In the model configuration section, it is mentioned that the model 3D structure is generated by perturbing a 1D background model from -10% to 10%. Given a 1D background model, how are number of layers decided per model example and how the perturbation and parameterization (using spherical harmonics) work to create a synthetic example?

2.	The results section is not very clear to understand especially due to lack of enough evidence for certain claims. I would highly appreciate if the paper could provide following information
(a)	Evidence to back claims in the results section related to training of H-Fourier models, how was physical constraint incorporated in the training, improve overall flow of forward and inverse results discussion with more figures/visualizations and share modeling related information
(b)	Number of examples generated in each tier
(c)	Incorporate more evaluation metrics such as MAE, MSE, SSIM. 

3.	The results section is primarily limited to acoustic dataset (tier-1) and provides very little insights on experimentation on the elastic dataset (tier-2). There is no discussion on training and evaluation results for the tier-3 dataset (real earth) which seems to be the main Earth’s scale dataset for global tomography. Would it be possible to have some visualizations of some datasets from all tiers along with experimental results for forward and inverse modeling for tier-2 and tier-3 dataset? 

4.	The paper explores model generalizability only for in-distribution cases, and it may be better to also study out-of-distribution generalizability to investigate ML effectiveness for both forward and inverse problems.

---

confidence: 5
contribution: 3
initial_rating: 6
presentation: 2
questions: See above
soundness: 3
strengths: The main strength of this work is in the data. There are many researchers in the scientific machine learning community that are looking for complex data sets that they can test their methods. I believe that the data set may be used for that.
summary: The paper propose a dataset for the global tomography problem. The data set consists of earth models and the wave fields that are obtained by using numerical simulations. The goal of this data is to be able to Tain a forward an an inverse network that allows for fast solutions of the global tomography problem.
weaknesses: There are a number of concerns that I have when considering this work.
1. The main problem is the diversity of the earth models. It is not clear to me that the models used are diverse enough to serve as a data set for the task at hand. The authors did not discuss this point clearly but this may be the most important point in the data set. I the models are not diverse enough then how can we learn something meaningful?

2. The second problem with the models is that we do not have ground truth and we will never have. Unlike problems in image processing, It is impossible to "open" the earth and obtain its physical properties at the resolution needed. How do you know you do not have bias in your models? 

3. The testing of the forward problem is not convincing. Wave equation solvers typically being judged on accuracy and on dispersion. Numerical dispersion is particularly important for solving inverse problem as it can introduce artifacts. It would be nice if the authors do not treat the wavefield as an image but rather follow most numerical PDE books when estimating the solution of a PDE.

4. How would the forward model do for models that are out of distribution? 

5. The inversion strategies use both optimization and direct methods where direct methods seems to win. Does this still work for out of distribution problems? 

6. What is the data fit for all the inversion methods. Do you actually fit the data? How does this compared with standard Gauss-Newton inversion where the data fit is being considered with the forward problem

Meta-review: The authors introduce the GlobalTomo dataset, which is a large 3D synthetic dataset tailored for seismic wavefield modeling and full-waveform tomography. For this problem there already exist a handful of publicly available datasets (such as OpenFWI and OpenFWI 2.0), The main novelty appears to be that while the previous datasets have focused on subsurface reconstruction and imaging, the proposed dataset includes "global-scale" wavefields (capturing core-to-surface phenomena at the planet scale). 

There is considerable computational effort involved in generating this dataset, but the authors do a good job of describing all the necessary details. Overall, the paper is well-written, and reviewers pointed out that the dataset can be valuable to the community. 

However, some concerns remained, particularly related to bias and generalization. There is no ground truth for the earth's interior; all existing models are inaccurate, and bias is inevitable in synthetically generated datasets; and performance may not generalize. Moreover, benchmark numbers on existing approaches may favor certain types of approaches over other ones. For these reasons I think the numerical results (and the rank ordering of the different inversion methods) reported in this manuscript need to be taken with appropriate caveats. 

This was a tough decision, but in the end the concern of generalizability is serious enough (and certainly not out of scope, as the authors argued in their response). The issue of sensitivity to OOD models warrants a more thorough discussion in the main paper with additional numerical results, over and above the checkerboard tests reported in the appendix. For maximizing impact, please consider using the best practices already established in the FWI community while constructing these additional studies.

Decision: Reject