Computational resources have increased in capacity over time — mostly by speed, partly by memory. Consequently, people have continuously explored the possibilities of performing wave modeling and inversion of increasing physical complexity. Achieving a detailed as possible image of the earth’s subsurface improves the success of hydrocarbon exploration, and it is important for other applications, such as archeology, mining, and engineering. I have developed an accurate computational method for elastic wave modeling up to tilted orthorhombic symmetry of anisotropy. The model may be covered by an arbitrary topographic function along the free surface. Through snapshots and seismograms of the wavefield, I confirm known effects from applying the code to plane, free surfaces (horizontal or tilted) as well as more complex topographies. The method is based on adapting a curved grid to a free-surface topography at hand, and transforming the wave equations and the topography free-surface boundary conditions from this grid to a rectangular grid, where finite-difference (FD) calculations can be performed. Free-surface topography boundary conditions for the particle velocities originate from locally setting the normal stress components to zero at the curved grid free surface. Vanishing normal traction is achieved by additionally imposing mirror conditions on stresses across the free surface. This leads me to achieve a more accurate modeling of free-surface waves (Rayleigh — Rg-waves in particular), using either FDs or any other numerical discretization method. Statics correction, muting, and destructive processing, which all consider free-surface effects as noise, can hence be avoided in inversion/imaging because surface effects can be more accurately simulated. By including near-surface effects in the full wavefield, we ultimately obtain superior inversion for interior earth materials, also for deeper physical medium properties.