Optimizing the finite-difference implementation of three-dimensional free-surface boundary in frequency-domain modeling of elastic waves

The problem of modeling seismic wave propagation for multiple sources, such as in the solution of gradient-based elastic full-waveform inversion, is an important topic in seismic exploration. The frequency-domain finite-difference (FD) method is a good choice for this purpose, mainly because of its simple discretization and high computational efficiency. However, when it comes to modeling the complete elastic wavefields, this approach has limited surface-wave accuracy because, when modeling with the strong form of the wave equation, it is not always easy to implement an accurate stress-free boundary condition. Although a denser spatial sampling is helpful for overcoming this problem, the additional discrete points will significantly increase the computational cost in the resolution of its resulting discrete system, especially in 3D problems. Furthermore, sometimes, when modeling with optimized schemes, an inconsistency in the computation precision between the regions at the free surface and inside the model volume would happen and introduce numerical artifacts. To overcome these issues, we have considered optimizing the FD implementation of the free-surface boundary. In our method, the problem was formulated in terms of a novel system of partial differential equations satisfied at the free surface, and the weighted-averaging strategy was introduced to optimize its discretization. With this approach, we can impose FD schemes for the free surface and internal region consistently and improve their discretization precision simultaneously. Benchmark tests for Lamb's problem indicate that the proposed free-surface implementation contributes to improving the simulation accuracy on surface waves, without increasing the number of grid points per wavelength. This reveals the potential of developing optimized schemes in the free-surface implementation. In particular, through the successful introduction of weighting coefficients, this free-surface FD implementation enables adaptation to the variation of Poisson's ratio, which is very useful for modeling in heterogeneous near-surface weathered zones.