Abstract Compressional (P) and shear (S) wave diffraction by free-surface topography plays a prominent part in the prediction of site responses for seismic risk estimation. Wave propagation modeling in 3-D media is required for an accurate estimation of these diffractions. We have extended the discrete wavenumber-indirect boundary integral equation method for a 3-D geometry in the case of irregular topography. The Green’s functions are expressed as finite sums of analytical density functions over the horizontal wavenumbers using the spatial periodicity of the topography and a discretization of the surface. We show that the evaluation over vertical wavenumber k z of the analytical integral is possible because a new factor in 1/ £ 2 exists. When the force point and the receiver point are at the same vertical position, we develop a numerical strategy to choose the sign of the exponential factor, which is not given by the analytical formulation. The free-streSs boundary conditions at the topography lead to a large linear system that can be solved to obtain the source density functions. Knowing these source density functions, we can compute the diffracted wave-field anywhere inside the medium. We have determined a useful optimal, imaginary frequency to obtain the displacement directly in the frequency domain, avoiding the necessity of returning to the time domain. We have then applied this method to investigate the effect of topography on the ground motion produced by a vertical incident P- or S-wavefield. Waveforms obtained for various topographic steepnesses and shapes show deterministic correlations between the maximum amplitude zone, the geometry of the 3-D topography, and the P- or 5-wave incident field characteristics. The maximum amplitude of the diffracted displacement is found near topographic zones that have horizontal or vertical dimensions closely related to the wavelength of the incident field. The predicted ground motion maximal amplifications are twice those calculated in the case of a flat topography.