Abstract We analyze the problem of a heterogeneous formulation of the equation of motion and propose a new 3D fourth-order staggered-grid finite-difference (FD) scheme for modeling seismic motion and seismic-wave propagation. We first consider a 1D problem for a welded planar interface of two half-spaces. A simple physical model of the contact of two media and mathematical considerations are shown to give an averaged medium representing the contact of two media. An exact heterogeneous formulation of the equation of motion is a basis for constructing the corresponding heterogeneous FD scheme. In a much more complicated 3D problem we analyze a planar-interface contact of two isotropic media (both with interface parallel to a coordinate plane and interface in general position in the Cartesian coordinate system) and a nonplanar-interface contact of two isotropic media. Because in the latter case 21 elastic coefficients at each point are necessary to describe the averaged medium, we consider simplified boundary conditions for which the averaged medium can be described by only two elastic coefficients. Based on the simplified approach we construct the explicit heterogeneous 3D fourth-order displacement-stress FD scheme on a staggered grid with the volume harmonic averaging of the shear modulus in grid positions of the stress-tensor components, volume harmonic averaging of the bulk modulus in grid positions of the normal stress-tensor components, and volume arithmetic averaging of density in grid positions of the displacement components. Our displacement-stress FD scheme can be easily modified into the velocity-stress or displacement-velocity-stress FD schemes. The scheme allows for an arbitrary position of the material discontinuity in the spatial grid. Numerical tests for 12 configurations in four types of models show that our scheme is more accurate than the staggered-grid schemes used so far. Numerical examples also show that differences in thickness of a soft surface or interior layer smaller than one grid spacing can cause considerable changes in seismic motion. The results thus underline the importance of having a FD scheme with sufficient sensitivity to heterogeneity of the medium.