Full-waveform inversion is a widely used method in seismic imaging, which involves the numerical modeling of the seismic wavefield. The question of how accurate this modeling is, is seldom discussed. This is a particular problem in land seismics and in the presence of topographic variations. For various reasons, the standard approach to waveform modeling in exploration seismology is the finite difference (FD) method. The implementation of a free surface boundary for a topographic surface in FD waveform modeling may be done as an internal boundary of the model in the heterogeneous approach or as an explicit condition at the free surface in the homogeneous approach. The widely used implementation of the image method (homogeneous), and the improved vacuum formulation (heterogeneous), are tested for their accuracy using two 2D FD modeling algorithms, ve2d_ref and DEINSE Black Edition, which implement these conditions to define an irregular free surface. These algorithms are first tested for their accuracy using analytical solutions in an elastic homogeneous full space, then for elastic homogeneous half-spaces with a planar and tilted free surface between 0 and 20°. Furthermore, the application to a viscoelastic inhomogeneous model with a topographic free surface is investigated. For all models, errors are introduced in the near-field for various reasons. At higher offsets, errors associated with the absorbing boundary frames or free surface become more prevalent. For a planar free surface, the accuracy tests show that the image method requires fewer points per minimum wavelength to accurately model the full waveforms than the improved vacuum formulation. In the case of a free surface not aligned with the axis of the FD grid, the free surface conditions show different behaviour, with the improved vacuum formulation producing higher errors than the image method at very small tilt angles (<5°). The accuracy of the improved vacuum formulation is highly dependent on the tilt angle. The convergence behaviour of the image method moves smoothly towards smaller errors, while the improved vacuum formulation remains constant and, after reaching a critical point in the discretisation, rapidly converges to small errors comparable, or better, than those achieved by the image method. This smooth convergence behaviour of the image method is easier to interpret than the fluctuation of relative differences in convergence tests, especially when dealing with inhomogeneous models with a topographic free surface, where there is no analytical solution available.