Finite-difference methods are popular for wave simulation within the seismic exploration community, thanks to their efficiency. However, difficulties arise when encountering complex topography due to the regular grid pattern of the finite-difference schemes. Despite alternatives that can handle the free surface with little effort, such as the spectral element or discontinuous Galerkin’s methods, incorporating a free-surface boundary condition within the finite-difference framework is still appealing, even at the cost of extra algorithm complexity and stronger requirement of computational resources. We present a free-surface boundary treatment within the finite-difference framework, belonging to the family of the immersed-boundary methods. Inherently, the presented boundary treatment is separated from the rest of the wave simulation, which makes it easy to be integrated in existing finite-difference codes. Specifically, we construct an extrapolation operator for each grid point above the free surface, if requested by the finite-difference stencil, to estimate its fictitious wavefield value at each time step. These operators are constructed only once and remain unchanged for all the time steps and source locations. The memory requirement of these operators is significant. Fortunately, grouping together multiple simulations concerning different source locations makes it possible to dilute the memory burden to a negligible level. Additionally, applying these operators incurs numerical noise, which may lead to long time instabilities. In such a scenario, additional numerical procedures, for instance, introducing artificial diffusion, are necessary to control the instability and obtain sensible simulation results. Successful applications of the presented boundary treatment to elastic-wave equations on domains with nontrivial topographies, in 2D and 3D, are presented. Robust and efficient numerical techniques to control high-frequency numerical noise remain to be investigated.