We examine a staggered pseudospectral method to solve a two-dimensional wav
e propagation problem with arbitrary nonlinear constitutive equations, and
evaluate a general image method to simulate the traction-free boundary cond
ition at the surface. This implementation employs a stress-velocity formula
tion and satisfies the free surface condition by explicitly setting surface
shear stress to zero and making the normal stress antisymmetric about the
free surface. Satisfactory agreement with analytical solutions to Lamb's pr
oblem is achieved for both vertical point force and explosion sources, and
with perturbation solutions for nonlinearly elastic wave propagation within
the domain of validity of such solutions. The Rayleigh wave, however, suff
ers much more severe numerical dispersion than do body waves. At four grids
per wavelength, the relative error in the Rayleigh-wave phase velocity is
25 times greater than the corresponding error in the body-wave phase veloci
ty. Thus for the Rayleigh wave, the pseudospectral method performs no bette
r than a low-order finite difference method.
A substantial merit of the image approach is that it does not assume any pa
rticular rheology, the method is readily applicable even when stresses are
not analytically related to kinematic variables, as is the case for most no
nlinear models. We use this scheme to investigate the response of a nonline
ar half-space with endochronic rheology, which has been fit to quasi-static
and dynamic observations. We find that harmonics of a monochromatic source
are generated and evolve with epicentral range, and energy is transferred
from low to higher frequencies for a broadband source. This energy redistri
bution characteristic of the propagation is strain-amplitude dependent, con
sistent with laboratory experiments. Compared with the linear response, the
nonlinear response of an endochronic layer near the surface shows a deampl
ification effect in the intermediate-frequency band and an amplification ef
fect in the higher-frequency band. The computational method, with modificat
ions to accommodate realistic nonlinear soil characteristics, could be appl
ied to estimate earthquake strong ground motions and path effects.