The 3-D numerical Chebychev modeling scheme accounts for surface topography. The method is based on spectral derivative operators. Spatial differencing in horizontal directions is performed by the Fourier method, whereas vertical derivatives are carried out by a Chebychev method that allows for the incorporation of boundary conditions into the numerical scheme. The method is based on the velocity‐stress formulation. The implementation of surface topography is done by mapping a rectangular grid onto a curved grid. Boundary conditions are applied by means of characteristic variables. The study of surface effects of seismic wave propagation in the presence of surface topography is important, since nonray effects such as diffractions and scattering at rough surfaces must be considered. Several examples show this. The 3-D modeling alogrithm can serve as a tool for understanding these phenomena since it computes the full wavefield.