A flexible and simple way of introducing stress-free boundary conditions for including three-dimensional (3D) topography in the finite-difference method is presented. The 3D topography is discretized in a staircase by stacking unit material cells in a staggered-grid scheme. The shear stresses are distributed on the 12 edges of the unit material cell so that only shear stresses appear on the free surface and normal stresses always remain embedded within the solid region. This configuration makes it possible to implement stress-free boundary conditions at the free surface by setting the Lamé coefficients λ and μ to zero without generating any physically unjustified condition. Arbitrary 3D topographies are realized by changing the distribution of λ and μ in the computational domain. Our method uses a parsimonious staggered-grid scheme that requires only 3/4 of the memory used in the conventional staggered-grid scheme in which six stress components and three velocity components need to be stored. Numerical tests indicate that 25 grids per wavelength are required for stable calculation. The finite-difference results are compared with those of the boundary-element method for the two-dimensional (2D) semi-circular canyon model. We also present the responses of a segment of semi-circular canyon and hemispherical cavity to vertically incident plane P, SV, and SH waves and discuss the response of a Gaussian hill to an isotropic point source embedded in the hill. In the segment of semi-circular canyon, the later portions of the synthetics are characterized by phases scattered from the two vertical side walls. The hemispherical cavity and 2D semi-circular canyon both show focusing of energy at the bottom of the cavity, although the focusing effect is stronger in the former geometry. Focusing and defocusing effects due to the strong topography of the Gaussian hill produce a strong amplification of displacements at a spot located on the flank opposite to the source. Backscattering from the top of the hill is also clearly seen.