Strongly imposing the free surface boundary condition for wave equations with finite difference operators

Acoustic and elastic wave equations are routinely used in geophysical and engineering studies to simulate the propagation of waves, with a broad range of applications, including seismology, near surface characterization, non-destructive structural evaluation, etc. Finite difference methods remain popular choices for these simulations due to their simplicity and efficiency. In particular, the family of finite difference methods based on the summation-by-parts operators and the simultaneous-approximation-terms technique have been proposed for these simulations, which offers great flexibility in addressing boundary and interface conditions. For the applications mentioned above, surface of the earth is usually associated with the free surface boundary condition. In this study, we demonstrate that the weakly imposed free surface boundary condition through the simultaneous-approximation-terms technique can have issue when the source terms, which introduces abrupt disturbances to the wave field, are placed too close to the surface. In response, we propose to build the free surface boundary condition into the summation-by-parts finite difference operators and hence strongly and automatically impose the free surface boundary condition to address this issue. The procedure is very simple for acoustic wave equation, requiring resetting a few rows and columns in the existing difference operators only. For the elastic wave equation, the procedure is more involved and requires special design of the grid layout and summation-by-parts operators that satisfy additional requirements, as revealed by the discrete energy analysis. In both cases, the energy conserving property is preserved. Numerical examples are presented to demonstrate the effectiveness of the proposed approach.