We investigate the parallel performance of an iterative solver for 3D heterogeneous Helmholtz problems related to applications in seismic wave propagation. For large 3D problems, the computation is no longer feasible on a single processor, and the memory requirements increase rapidly. Therefore, parallelization of the solver is needed. We employ a complex shifted-Laplace preconditioner combined with the Bi-CGSTAB iterative method and use a multigrid method to approximate the inverse of the resulting preconditioning operator. A 3D multigrid method with 2D semi-coarsening is employed. We show numerical results for large problems arising in geophysical applications.

The time-harmonic wave equation, also known as the Helmholtz equation, is obtained if the constant-density acoustic wave equation is transformed from the time domain to the frequency domain. Its discretization results in a large, sparse, linear system of equations. In two dimensions, this system can be solved efficiently by a direct method. In three dimensions, direct methods cannot be used for problems of practical sizes because the computational time and the amount of memory required become too large. Iterative methods are an alternative. These methods are often based on a conjugate gradient iterative scheme with a preconditioner that accelerates its convergence. The iterative solution of the time-harmonic wave equation has long been a notoriously difficult problem in numerical analysis. Recently, a new preconditioner based on a strongly damped wave equation has heralded a breakthrough. The solution of the linear system associated with the preconditioner is approximated by another iterative method, the multigrid method. The multigrid method fails for the original wave equation but performs well on the damped version. The performance of the new iterative solver is investigated on a number of 2D test problems. The results suggest that the number of required iterations increases linearly with frequency, even for a strongly heterogeneous model where earlier iterative schemes fail to converge. Complexity analysis shows that the new iterative solver is still slower than a time-domain solver to generate a full time series. We compare the time-domain numeric results obtained using the new iterative solver with those using the direct solver and conclude that they agree very well quantitatively. The new iterative solver can be applied straightforwardly to 3D problems.

SUMMARY Seismic surface wave analysis has provided important insight in the Earth's crustal and upper mantle structure and has recently become a standard tool in geotechnical engineering. Most current surface wave inversion methods are aimed at recovering near‐surface (shear) velocity profiles from dispersion curves, assuming a (smoothly varying) horizontally layered Earth. In some cases, however, one is interested in the location, strength or shape of local heterogeneities in the shallow subsurface. In this paper we focus on estimating the strength of near‐surface heterogeneity from scattered surface waves. This is a non‐linear inversion problem as the wavefield in the scatterer also depends on the contrast. For this reason the inversion is cast as an optimization problem in which we minimize the difference between the observed data and the modelled scattered field. The minimization problem is solved using a conjugate gradient algorithm. To accurately estimate the contrast we account for non‐linear interactions such as multiple scattering within the scattering domain in the forward problem. To do so, we use a domain‐type integral representation to express the near‐surface scattered wavefield, which is solved using the method of moments. The entire inversion is carried out in the frequency domain. This way, we may take advantage of the fact that the decay of surface waves with depth depends on frequency. We study 3‐D sensitivity kernels for the given inversion problem and observe that sensitivity of the wavefield with respect to heterogeneity in depth depends on frequency in the same way. We therefore, expect to be able to constrain heterogeneity in depth. A numerical example illustrates this and we conclude that it is in principle possible to reliably resolve heterogeneities and estimate their strengths. We compare the results from our algorithm to a similar but more efficient inversion scheme based on the Born approximation and show that, for the shallowest heterogeneities, this inversion can also recover the contrast well.

S U M M A R YWith the aid of a domain-type integral-equation formulation, the scattering problem is studied of an inclusion, embedded in a layered, traction-free, elastic half-space. This problem is characteristic for near-surface scattering problems occurring in exploration seismology. The singular character of the Green's function is taken into account in an accurate and efficient manner using an asymptotic description of its near-field behaviour. The numerical results show good qualitative agreement with experimental scale-model data. From a number of model studies, it is concluded that we can accurately model near-surface scattering effects, which can seriously distort the wave fronts of upcoming reflections.