We examine a staggered pseudospectral method to solve a two-dimensional wave propagation problem with arbitrary nonlinear constitutive equations, and evaluate a general image method to simulate the traction-free boundary condition at the surface. This implementation employs a stress-velocity formulation 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 problem 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, suffers 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 velocity. Thus for the Rayleigh wave, the pseudospectral method performs no better than a low-order finite difference method.
A substantial merit of the image approach is that it does not assume any particular rheology, the method is readily applicable even when stresses are not analytically related to kinematic variables, as is the case for most nonlinear models. We use this scheme to investigate the response of a nonlinear 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 redistribution characteristic of the propagation is strain-amplitude dependent, consistent with laboratory experiments. Compared with the linear response, the nonlinear response of an endochronic layer near the surface shows a deamplification effect in the intermediate-frequency band and an amplification effect in the higher-frequency band. The computational method, with modifications to accommodate realistic nonlinear soil characteristics, could be applied to estimate earthquake strong ground motions and path effects.