MR elastography (MRE) is an MRI modality that is increasingly being used to image tissue elasticity throughout the body. One MRE technique that has received a great deal of attention is based on visualizing shear waves, which reveal stiffness by virtue of their local wavelength. However, the shape of propagating shear waves can also provide valuable information about the nonlinear stress-strain behavior of tissue. Here an experiment is proposed that allows the observation of nonlinear wave propagation based on spatial-temporal phase contrast images. A theoretical description of the wave propagation was developed that reflects typical MRE excitation, which involves excitation modes both parallel and perpendicular to B0. Based on this model, it is shown that both odd and even higher harmonics are produced with their amplitudes dependent on the details of the actuator, imaging geometry, and the nonlinear tissue properties. With appropriate motion encoding, harmonic vibrations arising from nonlinear tissue response can be detected. The effect is demonstrated on an agarose gel phantom using a sinusoidal shear vibration of 150 Hz, and clearly shows the presence of harmonics at 600 and 750 Hz. Using an estimate of the strain energy of the phantom, we were able to determine the nonlinear tissue properties.