An accurate free-surface boundary condition is important for solving a wide variety of seismic modeling problems. In particular, for earthquake site studies or shallow environmental investigations the surface of the earth may have a significant impact on the outcome of simulations. Computations based on several elastic/viscoelastic flat horizontal free-surface conditions are compared and benchmarked against an analytical solution. An accurate and simple condition is found and then generalized to allow for irregular free surfaces. This new method is simple to implement in conventional staggered finite-difference schemes, is computationally efficient and enables modeling of highly irregular topography. The accuracy of the method is investigated and criteria for sampling of the wavefield are derived.