Numerical simulation of P-wave propagation in 2D VTI media with a topographic free surface by a curvilinear-grid finite-difference method

Topographic free surfaces and seismic anisotropy significantly modify seismic wave propagation. With increasing demands for high-precision subsurface imaging, advanced seismic prospecting and seismological imaging methods require simultaneous consideration of topographic free surface and anisotropy effects on P-wave propagation. However, conventional finite-difference methods struggle to accurately simulate both topographic free surface and anisotropy in first-order acoustic velocity-stress equations, failing to capture complex wavefield responses. To address this challenge, we extend the existing curvilinear-grid finite-difference method, originally developed for the elastic wave VTI system of velocity-stress equations, to the acoustic VTI velocity-stress system with topographic free surfaces. The primary difficulty lies in setting the shear-wave velocity to zero in the elastic wave equations, which causes the free-surface velocity-derivative constraint matrix for acoustic topographic free-surface models to become rank-deficient or severely ill-conditioned. Therefore, we derive the velocity-derivative constraint specific to the acoustic VTI system and use the Moore-Penrose generalized inverse to obtain the normal velocity derivatives needed for updating the free-surface boundary condition. The resulting formulation includes body-fitted curvilinear grids, collocated-grid DRP/opt MacCormack discretization, and traction-image boundary treatment. Numerical experiments on flat-surface, Gaussian-topography, and BP benchmark models demonstrate stable performance under the tested conditions, while near free-surface waveforms computed by our approach closely match independent reference solutions. This method provides a reliable forward-modeling foundation for 2D acoustic-wave applications that require simultaneous, high-accuracy consideration of irregular topography and anisotropy.