  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 U-shaped
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.
