Abstract This article provides a technique to model seismic motions in 3D elastic media using fourth-order staggered-grid finite-difference (FD) operators implemented on a mesh with nonuniform grid spacing. The accuracy of the proposed technique has been tested through comparisons with analytical solutions, conventional 3D staggered-grid FD with uniform grid spacing, and reflectivity methods for a variety of velocity models. Numerical tests with nonuniform grids suggest that the method allows sufficiently accurate modeling when the grid sampling rate is at least 6 grid points per shortest shear wavelength. The applicability for a finite fault with nonuniform distribution of point sources is also confirmed. The use of nonuniform spacing improves the efficiency of the FD methods when applied to large-scale structures by partially avoiding the spatial oversampling introduced by the uniform spacing in zones with high velocity. The significant reduction in computer memory that can be obtained by the new technique improves the efficiency of the 3D-FD method at handling shorter wavelengths, larger areas, or more realistic 3D velocity structures.