We have implemented and verified a parallel‐series Iwan‐type nonlinear model in a 3D fourth‐order staggered‐grid velocity–stress finite‐difference method. The Masing unloading and reloading behavior is simulated by tracking an overlay of concentric von Mises yield surfaces. Lamé parameters and failure stresses pertaining to each surface are calibrated to reproduce the stress–strain backbone curve, which is controlled by the reference strain assigned to a given depth level. The implementation is successfully verified against established codes for 1D and 2D SH‐wave benchmarks. The capabilities of the method for large‐scale nonlinear earthquake modeling are demonstrated for an 7.8 dynamic rupture ShakeOut scenario on the southern San Andreas fault. Although ShakeOut simulations with a single yield surface reduces long‐period ground‐motion amplitudes by about 25% inside a waveguide in greater Los Angeles, Iwan nonlinearity further reduces the values by a factor of 2. For example, inside the Whittier Narrows corridor spectral accelerations at a period of 3 s are reduced from 1g in the linear case to about 0.8 in the bilinear case and to 0.3–0.4g in the multisurface Iwan nonlinear case, depending on the choice of reference strain. Normalized shear modulus reductions reach values of up to 50% in the waveguide and up to 75% in the San Bernardino basin at the San Andreas fault. We expect the implementation to be a valuable tool for future nonlinear 3D dynamic rupture and ground‐motion simulations in models with coupled source, path, and site effects.