Full Waveform Inversion Research Articles

Review of physics-informed machine-learning inversion of geophysical data

We review five types of physics-informed machine-learning (PIML) algorithms for inversion and modeling of geophysical data. Such algorithms use the combination of a data-driven machine-learning (ML) method and the equations of physics to model or invert geophysical data (or both). By incorporating the constraints of physics, PIML algorithms can effectively reduce the size of the solution space for ML models, enabling them to be trained on smaller data sets. This is especially advantageous in scenarios in which data availability may be limited or expensive to obtain. In this review, we restrict the physics to be that from the governing wave equation, either as a constraint that must be satisfied or by using numerical solutions of the wave equation for modeling and inversion. This approach ensures that the resulting models adhere to physical principles while leveraging the power of ML to analyze and interpret complex geophysical data. There are several potential benefits of PIML compared to standard numerical modeling or inversion of seismic data computed by, for example, finite-difference solutions to the wave equation. 1) Empirical tests suggest that PIML algorithms constrained by the physics of wave propagation can sometimes resist getting stuck in a local minima compared with standard full-waveform inversion (FWI). 2) After the weights of the neural network are found by training, then the forward and inverse operations by PIML can be more than several orders of magnitude more efficient than FWI. However, the computational cost for general training can be enormous. 3) If the ML inversion operator [Formula: see text] is locally trained on a small portion of the recorded data [Formula: see text], then there is sometimes no need for millions of training examples that aim for global generalization of [Formula: see text]. The benefit is that the locally trained [Formula: see text] can be used to economically invert the remaining test data [Formula: see text] for the true velocity [Formula: see text], where [Formula: see text] can comprise more than 90% of the recorded data.

WUS324: Multiscale Full Waveform Inversion Approaching Convergence Improves Waveform Fits While Imaging Seismic Structure of the Western United States

AbstractWe report a new model of radially anisotropic crustal and upper mantle structure of the western United States (WUS324) obtained from full waveform inversion of earthquake data. We ran three multiscale inversion stages beyond model WUS256 (Rodgers et al., 2022, https://doi.org/10.1029/2022jb024549) allowing them to approach convergence to fit a larger data set to a shorter minimum period of 16 s. WUS324 is based on 324 total iterations from its starting model, significantly more (16 times) than previous studies. Waveform misfit reductions are 66%–70% for both the inversion data and an independent validation data set providing confidence in the predictive power of the model. WUS324 provides much better fits and reveals shear wavespeed, vS, structure of this large region with more detail than previous waveform tomography models. We show representative images demonstrating the resolution of diverse seismic structure across this highly heterogeneous region including oceanic lithosphere, subducting slabs and continental magmatism.

Physics-guided full waveform inversion using Encoder-Solver convolutional neural networks

Abstract Full Waveform Inversion (FWI) is an inverse problem for estimating the wave velocity distribution in a given domain, based on observed data on the boundaries. The inversion is computationally demanding because we are required to solve multiple forward problems, either in time or frequency domains, to simulate data that are then iteratively fitted to the observed data. We consider FWI in the frequency domain, where the Helmholtz equation is used as a forward model, and its repeated solution is the main computational bottleneck of the inversion process. To ease this cost, we integrate a learning process of an Encoder-Solver preconditioner that is based on convolutional neural networks (CNNs). The Encoder-Solver is trained to effectively precondition the discretized Helmholtz operator given velocity medium parameters. Then, by re-training the CNN between the iterations of the optimization process, the Encoder-Solver is adapted to the iteratively evolving velocity medium as part of the inversion. Without retraining, the performance of the solver deteriorates as the medium changes. Using our light retraining procedures, we obtain the forward simulations effectively throughout the process. We demonstrate our approach to solving FWI problems using 2D geophysical models with high-frequency data.

FWI-WF: Robust full-waveform inversion using hybrid Wasserstein-1 and Fourier metrics

The performance of full-waveform inversion (FWI) in constructing high-resolution subsurface models is closely related to the design of mismatch functions. The least-squares norm ([Formula: see text]) is commonly used; however, it is prone to local minima when high-quality initial guess and low-frequency data are unavailable. The Wasserstein-1 ([Formula: see text]) metric captures time shifts more effectively, but it may be plagued by imprecise deep structures. The Fourier metric leverages power spectra from simulated and observed data, offering higher-resolution updates near solutions. We develop a progressive waveform inversion method called FWI-WF using [Formula: see text] and Fourier metrics. Specifically, in the early stage of inversion, we apply greater weight to the [Formula: see text] metric for constructing a good background model and avoiding falling into local minima. Then, the Fourier metric gradually dominates to refine edges and deep structures, providing high-resolution inversion results. During the optimization process, we use automatic differentiation to improve inversion efficiency. Experimental results on three baseline geologic models indicate that FWI-WF outperforms three state-of-the-art methods.

Study on non-destructive visualization identification of concrete internal cracks based on SH array ultrasound

Ultrasonic nondestructive testing (NDT) technology was one of the most frequently used and fastest-developing NDT techniques in detection field. This method had high detection sensitivity, more accurate defect localization, and was harmless to the human body. With the development of ultrasonic testing technology, modern ultrasonic imaging technology generally had automatic data acquisition and processing functions, which allowed for intuitive and clear recognition of internal defects by the human eyes. The overall purpose of this study was to visualize the cracks in concrete by synthetic aperture focusing technique (SAFT) and full-waveform inversion (FWI) algorithms on the basis of ultrasonic nondestructive testing technology. At present, there have been studies on the use of ultrasonic instruments to detect the internal defects of concrete. However, the accuracy and authenticity of this method have not been compared in previous studies. Therefore, in this experiment, four concrete specimens with different buried depths, angles, shapes, and thicknesses of cracks were designed. Based on SH ultrasonic wave NDT technology, two types of operations were performed on the collected ultrasonic wave data: one was the application of SAFT method, and the other was the first application of FWI method for detecting internal cracks in concrete. The study found that both algorithms were most accurate in determining the burial depth of cracks, with error results not exceeding 10 %. They also showed good recognition effects for the angle, shape, and thickness of the buried cracks. This proved that both identification methods had the capability to recognize internal cracks. The Introview software bundled with the ultrasonic instrument could identify cracks using SAFT directly and quickly. When FWI was used with MATLAB computation, the results were more precise and intuitive. The overall purpose of this study was to detect cracks in concrete by SAFT and FWI algorithms on the basis of ultrasonic nondestructive testing technology. Combining the two crack identification techniques could be of great significance for ensuring the long-term safety and stability of concrete structures, providing strong support for structural maintenance and reinforcement.

A comprehensive crosshole seismic experiment in glacial sediments at the ICDP DOVE site in the Tannwald Basin

Abstract. Glaciers have shaped the Alpine landscape by carving deep valleys and depositing sediments to form overdeepened basins. Understanding these processes provides information on the evolution of the climate and landscape. One such overdeepened structure is the Tannwald Basin (ICDP site 5068_1) north of Lake Constance, which was formed by the Rhine Glacier in several glacial cycles. In order to study these sediments and their seismic properties down to about 160 m depth, we conducted seismic crosshole experiments between three boreholes, obtaining compressional (P) wave data. The P-wave data are generated by a sparker source and recorded by a 24-station hydrophone string. We present the data acquisition and review our approach for future optimization, suggesting a finer time sampling interval and a separate registration of borehole and surface receivers. Travel-time tomography of the P-wave first-arrival picks under geostatistical constraints yields initial subsurface models. The tomograms correlate well with cased-hole sonic logs and the lithology derived from the core of one of the boreholes. These results will be further investigated in future research, which will include full-waveform inversion (FWI) to obtain high-resolution subsurface models.

Full-Waveform Inversion Imaging of Cortical Bone Using Phased Array Tomography.

Classic ultrasound bone imaging modalities usually demand either a prior knowledge or an advanced estimation on speed of sound (SoS), which not only renders to a burdensome imaging process but also supplies a limited resolution. To overcome these drawbacks, this article proposed a frequency-domain full-waveform inversion (FDFWI) modality using phased array tomography for high-accuracy cortical bone imaging. A transmission scenario of ultrasound wave in 2-D space was presented in the frequency domain to simulate the forward wavefield propagation. Iterations in the inversion process were performed by matching the simulation wavefield to the experimental one from low to high discrete frequency points. Moreover, the association between the maximum initial frequency and the initial SoS model was explored to prevent the occurrence of cycle-skipping phenomenon, which could lead to the outcomes being trapped in local minima. The feasibility and effectiveness of the proposed imaging scheme were testified by simulation, phantom, and ex-vivo studies, with mean relative errors of cortical part being 3.18%, 8.71%, and 9.36%, respectively. It is verified that the proposed FDFWI method is an effective way for parametric imaging of cortical bone without any prior knowledge of sound speed.

Multi‐Method Structural Investigation of the Schneiderberg–Baalberge Burial Mound (Saxony‐Anhalt, Germany) Including Seismic Full‐Waveform Inversion (FWI)

ABSTRACTThe construction history and subsequent usage of burial mounds are an important testimony for socio‐economic transformation in prehistoric societies. The Baalberge–Schneiderberg burial mound, subject of the presented study, falls in this category as it is considered as an important monument that indicates the emergence of early social stratification during the Chalcolithic period in central Europe. This hypothesis relies on the chronological development of the burial mound, which is not fully understood until now. Therefore, a reconstruction of the complex stratigraphy of the burial mound including construction phases and later alterations is highly relevant for archaeological research, but the required excavations would be onerous and inconsistent with preservation efforts. In this paper, we demonstrate that non‐invasive geophysical prospection, especially seismic sounding with shear and Love waves, is suitable to obtain the required stratigraphic information, if seismic full waveform inversion (FWI) and reflection imaging are applied. Complementary information on the preservation state of the mound is obtained through Electrical Resistivity Tomography (ERT) and Electromagnetic Induction (EMI) measurements. To support the seismic and geoelectric results, we utilize Dynamic Testing (DynP), geoarchaeological corings, 14C‐Dating and archaeological records. Our investigations reveal two construction phases of the Baalberge–Schneiderberg mound. The 14C‐Dating yields dates for the older burial mound that are contemporary to the Chalcolithic Baalberge group (4000–3400 bc). During the Early Bronze Age (EBA), the mound was enlarged to its final size by people of the Aunjetitz/Únětice society (2300–1600 bc). However, both seismic and geoelectric depth sections show an extensive disturbance of the original stratigraphy due to former excavations. For this reason, the exact shape of the older burial mound cannot be determined exactly. Based on our data, we estimate that its height was below 2 m. In consequence, the original Baalberge burial mound was less monumental as until now assumed, which potentially prompting a revision of its significance as indicator for social differentiation.

Optimising seismic imaging design parameters via bilevel learning

Abstract Full waveform inversion (FWI) is a standard algorithm in seismic imaging. It solves the inverse problem of computing a model of the physical properties of the earth’s subsurface by minimising the misfit between actual measurements of scattered seismic waves and numerical predictions of these, with the latter obtained by solving the (forward) wave equation. The implementation of FWI requires the a priori choice of a number of ‘design parameters’, such as the positions of sensors for the actual measurements and one (or more) regularisation weights. In this paper we describe a novel algorithm for determining these design parameters automatically from a set of training images, using a (supervised) bilevel learning approach. In our algorithm, the upper level objective function measures the quality of the reconstructions of the training images, where the reconstructions are obtained by solving the lower level optimisation problem—in this case FWI. Our algorithm employs (variants of) the BFGS quasi-Newton method to perform the optimisation at each level, and thus requires the repeated solution of the forward problem—here taken to be the Helmholtz equation. This paper focuses on the implementation of the algorithm. The novel contributions are: (i) an adjoint-state method for the efficient computation of the upper-level gradient; (ii) a complexity analysis for the bilevel algorithm, which counts the number of Helmholtz solves needed and shows this number is independent of the number of design parameters optimised; (iii) an effective preconditioning strategy for iteratively solving the linear systems required at each step of the bilevel algorithm; (iv) a smoothed extraction process for point values of the discretised wavefield, necessary for ensuring a smooth upper level objective function. The algorithm also uses an extension to the bilevel setting of classical frequency-continuation strategies, helping avoid convergence to spurious stationary points. The advantage of our algorithm is demonstrated on a problem derived from the standard Marmousi test problem.

Seismic evidence of an extremely thin crust along a 70 Ma slow-spreading crustal segment in the equatorial Atlantic Ocean

The oceanic crust at slow-spreading ridges is formed by a combination of magmatic and tectonic processes. Seismic tomographic studies, commonly used to determine crustal structures, indicate that mature slow-spreading crust has an average thickness of 6.1 ± 0.9 km, with a lower to upper crustal thickness ratio of 2.3. However, tomographic methods are based on high-frequency approximations and can only provide information about large-scale crustal structures. Using seismic full-waveform inversion applied to ocean bottom seismometer data, we report the presence of nearly uniform and extremely thin crust (4 ± 0.2 km) with a lower to upper crustal thickness ratio of 0.74 along an ∼100 km, 70 Ma crustal segment formed at the slow-spreading Mid-Atlantic Ridge in the equatorial Atlantic Ocean. This thin crust is underlain by a Moho transition zone that is 1 ± 0.5 km thick, resulting in a combined crustal and Moho transition zone thickness of 5 ± 0.6 km. The crustal thinning is primarily due to the thinning of the lower crust, suggesting that a significant part of the crust is formed by lava flows and dike intrusions. The presence of an extremely thin crust indicates a decrease in mantle temperature by at least 50 °C, suggesting a vast extent of the upper mantle thermal minimum in the equatorial Atlantic Ocean, and the absence of influence of the mantle from the Siera Leone mantle plume 70 m.y. ago.

Assessing the value of combined use of distributed acoustic sensing (DAS) and multicomponent vertical seismic profiling (VSP) data with network-based full waveform inversion and uncertainty quantification: A case study in Alberta, Canada

Distributed acoustic sensing (DAS) technology, deployed in a vertical seismic profiling (VSP) experimental configuration, has emerged as a candidate for non-disruptive and low-cost seismic monitoring of CO2 geostorage and plume evolution. Full waveform inversion (FWI) has likewise received significant attention because it uses relatively complete physical models of wave propagation and because of its sample-by-sample incorporation of data information. Recent artificial neural network-based FWI algorithms (built with, for instance, recursive neural networks, or RNNs) have added to FWI a range of flexible and efficient tools for gradient computation and options for uncertainty assessment and initial model proxies. An important current research area for the use of DAS data is to better understand how they change our confidence levels in models selected through FWI. In particular, we seek to understand whether either DAS data or conventional geophone data alone are optimal for FWI model selection in the CO2 problem, and if not, to what degree they complement each other. The Snowflake 4D VSP dataset, which includes multi-offset and multi-azimuth broadband sources illuminating both fiberoptic cable and densely sampled accelerometer phones in the borehole, was acquired by our group to directly address these questions. In this study, we quantify uncertainty (UQ) by sampling the posterior model covariance matrix from the inverse Hessian matrix at the end of RNN-FWI runs on the Snowflake baseline data, invoking a velocity-density parameterization, and involving mixtures of accelerometer and DAS data. In this UQ context, the complementary effect of combining accelerometer and DAS data is evident in the Vp and Vs models. Our confidence in the approach is bolstered by its independent prediction of a region of known high uncertainty. In the overall pursuit of reliable and low-cost monitoring tools, this supports continued consideration of a multicomponent sensors supported by DAS approach.

3D structural modelling of the Kopili fault zone in North-East India: Seismotectonic analysis utilising focal mechanism solutions of small-to-moderate earthquakes

This work presents new findings from recent seismic data which have never been exploited so far for the weak-to-moderately energetic (ML 2.3–6.4) earthquakes and to depict a 3D scheme for the Kopili Fault (KF) Zone. For the first time, we figure out different fault kinematics in the upper crust such as reverse and strike-slip mechanisms located along the KF zone. Moment Tensor (MT) solutions are successfully obtained through full-waveform inversions for 60 seismic events. The stress inversions are estimated based on the MT results as well as previous findings from 16 small to moderate earthquakes (MW 3.6–6.3). The MT solutions demonstrate that strike-slip kinematics is the most prevalent while a few thrust/normal earthquake events are observed in the area. The seismic cross-sections display that these earthquakes took place at a depth of 40 km and are mostly generated on the eastern side of the fault. According to the 3D spatial and temporal distribution of relocated events, the KF zone is not a completely vertical fault but rather one with a slight eastward dip. Our stress inversion results also indicate that this fault is equivalent to a fault plane with a dip angle of around 80˚, which is compatible with the previously estimated dip angle (75˚). We interpret that the research region is under a transpressional stress regime, which is likely to be the reason for the current activity of the KF zone. We postulate that our findings can be useful to enhance seismic hazard assessment in the densely populated region, which is also covered by vital infrastructure.

A global approach to detecting and characterizing water leakage in a concrete bridge deck: Parametric study to validate an adapted full-waveform inversion method

The objective of this study is to address issues related to defects in waterproofing membranes through the use of Ground Penetration Radar to detect water leakage into the concrete layer of bridge decks. Given that processing radar signals based solely on temporal data introduces significant estimation errors, an advanced method optimized from Full-Waveform Inversion (FWI) is employed. This method accounts for several unknown factors, including boundary conditions, antenna positioning inaccuracies, and approximations inherent in 2D modeling during the inversion process. The method is adaptable and capable of reconstructing both the dielectric and geometric parameters of a multilayered structure with any number of layers and unknown parameters using just a few A-scans (up to 10, depending on the number of parameters and layers). In contrast, other methods typically require several hundred A-scans. This efficiency is achieved due to the two-dimensional nature of the layer system. Additionally, the simplicity of the structure facilitates a much faster and more straightforward inversion algorithm, hence making these features especially advantageous for practical applications. The optimized Full-Waveform Inversion approach allows for an accurate determination of unknown parameters within a multilayered medium. The high accuracy of this method is validated through a direct comparison with experiments, wherein the exact parameters are known. Such an approach enables the attainment of a low relative error using the currently available measurement device.

Seismic Evidence for Velocity Heterogeneity Along ∼40 Ma Old Oceanic Crustal Segment Formed at the Slow‐Spreading Equatorial Mid‐Atlantic Ridge From Full Waveform Inversion of Ocean Bottom Seismometer Data

AbstractIn slow spreading environments, oceanic crust is formed by a combination of magmatic and tectonic processes. Using full waveform inversion applied to active‐source ocean bottom seismometer data, we reveal the presence of a strong lateral variability in the 40–48 Ma old oceanic crust formed at the slow‐spreading Mid‐Atlantic Ridge in the equatorial Atlantic Ocean. Over a 120 km‐long section between the St Paul fracture zone (FZ) and the Romanche transform fault (TF), we observe four distinct 20–30 km long crustal segments. The segment affected by the St Paul FZ consists of three layers, an ∼2 km thick layer with a P‐wave velocity <6 km/s, a 1.5 km thick middle crust with a velocity of 6–6.5 km/s, and an underlying layer where velocity is ∼7 km/s, representing the lower crust. The segment associated with an abyssal hill morphology contains a high velocity of ∼7 km/s at 2–2.5 km below the basement, indicating the presence of primitive gabbro or serpenized peridotite. The segment associated with a low basement morphology seems to have 5.5–6.5 km/s velocity starting near the basement extending down to ∼4 km depth, indicating chemically distinct crust. The segment close to the Romanche TF, a velocity 4.5–5 km/s near the seafloor increasing to 7 km/s at 4 km depth indicates a magmatic origin. The four distinct crustal segments have a good correlation with the overlying seafloor morphology. These observed strong crustal heterogeneities could result from alternate tectonic and magmatic processes along the ridge axis, possibly modulated by thermal and/or chemical variations in the mantle during their formation along the ridge segment.

Modeling Rayleigh wave in viscoelastic media with constant Q model using fractional time derivatives

The propagation of seismic waves within the near-surface weathering layers, characterized by their low-quality factors (Q), is often accompanied by strong attenuation and dispersion phenomena. Among these, the Rayleigh wave, with its sensitivity to dispersion, has proven to be a powerful tool for near-surface exploration. We propose a novel approach for simulating Rayleigh wave propagation in such low-Q media. Our method uses the time-domain fractional wave equation with memory effect, based on Kjartansson's constant-Q (CQ) model, for accurate characterization of the propagation process. To solve numerically the wave equation with the fractional derivatives, we employ a finite-difference method combined with the auxiliary differential equation-perfectly matched layer (ADE-PML) and the acoustic-elastic boundary approach (AEA). The algorithm's high computational accuracy is verified through comparison with the conventional integer-order wave equation based on the nearly constant-Q (NCQ) models in strong attenuation media. The research in this paper deepens our understanding of the propagation characteristics of Rayleigh waves in strongly weathering layers. This new method strongly supports those seismic imaging and inversion methods depending on seismic modeling, including the reverse time migration and the full waveform inversion of the internal structure of low-Q media.

Enhancing Full Waveform Inversion of Field GPR Data: A Source-Independent Approach with Dynamic Reference Selection via SE-Wave-U-Net

Ground Penetrating Radar (GPR) is a widely-utilized non-destructive detection tool. Full Waveform Inversion (FWI) offers a reliable algorithm for accurately imaging GPR data. Despite its potential, FWI's effectiveness is often compromised by unknown excitation sources, which can lead to computational failures and/or diminished imaging accuracy. Addressing this challenge, this paper introduces a convolution-type objective function designed to mitigate the detrimental effects of these unknown excitation sources on FWI imaging. Critical to the success of this function is the precise selection of reference traces. Further refining the capability for signal separation, we have integrated Squeeze-and-Excitation (SE) modules into the traditional Wave-U-Net architecture, culminating in the development of the SE-Wave-U-Net framework. This framework is adept at dynamically extracting accurate reference traces from field GPR data, which are essential for the accurate FWI imaging. The efficacy, accuracy, and practical applicability of proposed methods are validated through extensive analysis of numerical, laboratory experiments, and field GPR data. The proposed algorithm not only establishes a high-precision computational framework for FWI but also enhances the reliability of GPR imaging in complex environments.

Greedy selection of optimal location of sensors for uncertainty reduction in seismic moment tensor inversion

We address an optimal sensor placement problem through Bayesian experimental design for seismic full waveform inversion for the recovery of the associated moment tensor. The objective is that of optimally choosing the location of the sensors (stations) from which to collect the observed data. The Shannon expected information gain is used as the objective function to search for the optimal network of sensors. A closed form for such objective is available due to the linear structure of the forward problem, as well as the Gaussian modeling of the observational errors and prior distribution. The resulting problem being inherently combinatorial, a greedy algorithm is deployed to sequentially select the sensor locations that form the best network for learning the moment tensor. Numerical results are presented to display the optimal network of sensors and how the uncertainty in the inferred seismic moment tensor contracts. The scenario of full three-dimensional velocity models or unknown earthquake source locations is treated as nuisance uncertainty, contributing to the overall uncertainty without being the focus of the inversion. This is addressed using a consensus approach over a set of realizations of the nuisance parameter. We analyzed the resulting network of stations for the moment tensor inversion under model misspecification, which reflects realistic data-generating processes.

Full waveform inversion with smoothing of dilated convolutions

Full waveform inversion with smoothing of dilated convolutions

Adjoint-based inversion for stress and frictional parameters in earthquake modeling

We present an adjoint-based optimization method to invert for stress and frictional parameters used in earthquake modeling. The forward problem is linear elastodynamics with nonlinear rate-and-state frictional faults. The misfit functional quantifies the difference between simulated and measured particle displacements or velocities at receiver locations. The misfit may include windowing or filtering operators. We derive the corresponding adjoint problem, which is linear elasticity with linearized rate-and-state friction and, for forward problems involving fault normal stress changes, nonzero fault opening, with time-dependent coefficients derived from the forward solution. The gradient of the misfit is efficiently computed by convolving forward and adjoint variables on the fault. The method thus extends the framework of full-waveform inversion to include frictional faults with rate-and-state friction. In addition, we present a space-time dual-consistent discretization of a dynamic rupture problem with a rough fault in antiplane shear, using high-order accurate summation-by-parts finite differences in combination with explicit Runge–Kutta time integration. The dual consistency of the discretization ensures that the discrete adjoint-based gradient is the exact gradient of the discrete misfit functional as well as a consistent approximation of the continuous gradient. Our theoretical results are corroborated by inversions with synthetic data. We anticipate that adjoint-based inversion of seismic and/or geodetic data will be a powerful tool for studying earthquake source processes; it can also be used to interpret laboratory friction experiments.

Deep neural helmholtz operators for 3D elastic wave propagation and inversion

Summary Numerical simulations of seismic wave propagation in heterogeneous 3D media are central to investigating subsurface structures and understanding earthquake processes, yet are computationally expensive for large problems. This is particularly problematic for full waveform inversion, which typically involves numerous runs of the forward process. In machine learning there has been considerable recent work in the area of operator learning, with a new class of models called neural operators allowing for data-driven solutions to partial differential equations. Recent works in seismology have shown that when neural operators are adequately trained, they can significantly shorten the compute time for wave propagation. However, the memory required for the 3D time domain equations may be prohibitive. In this study, we show that these limitations can be overcome by solving the wave equations in the frequency domain, also known as the Helmholtz equations, since the solutions for a set of frequencies can be determined in parallel. The 3D Helmholtz neural operator is 40 times more memory-efficient than an equivalent time-domain version. We employ a Helmholtz neural operator for 2D and 3D elastic wave modeling, achieving two orders of magnitude acceleration compared to a baseline spectral element method. The neural operator accurately generalizes to variable velocity structures and can be evaluated on denser input meshes than used in the training simulations. We also show that when solving for wavefields strictly on the surface, the accuracy can be significantly improved via a graph neural operator layer. In leveraging automatic differentiation, the proposed method can serve as an alternative to the adjoint-state approach for 3D full-waveform inversion, reducing the computation time by a factor of 350.