Four-level two-dimensional (213) Fourier transform relaxation spectra are simulated with response functions for a chromophore pair in the exponential relaxation (optical Bloch model) limit. The parameters in this study are chosen to model coupled carbonyl stretching vibrations. As long as coherence persists, every peak in the real 2D spectra has a partially mixed absorptive/dispersive ("phase-twisted") shape because the nonlinear signals are not symmetric with respect to interchange of the first two pulses. This asymmetry in 2D relaxation spectra arises from coherence between singly excited states and a red shift of the doubly excited state. Coherence between the singly excited states causes oscillation of the 2D spectra and the associated spectrally resolved pump-probe (SRPP) transients at the quantum beat frequency. Projecting the phase-twisted nature of the 2D peaks onto the detection frequency axis, the SRPP peaks are also asymmetric about their maximum when not at maximum or minimum amplitude. Three-dimensional Fourier transform (3DFT) methods are used to simulate absorption/dispersion and beam geometry distortions of the multilevel 2D spectra with cross peaks. The distortions can be understood by consideration of their effects on individual coherence pathways that contribute to peaks in the 2D spectra. The beam geometry distortion explains some unequal cross peak amplitudes previously observed experimentally by Khalil et al. (J. Chem. Phys. 2004, 121, 362). A representation of 2D spectra that reduces beam geometry distortion is presented. If the transformation to correct for beam geometry distortion is combined with the transformations that correct absorptive/dispersive propagation distortions (J. Chem. Phys. 2007, 126, 044511), the recovered 2D spectrum matches the ideal 2D spectrum after all coherence is destroyed. In the presence of coherence, the new representation reduces the error in the distorted 2D spectrum by a factor of 4 for practical 2D-IR experimental conditions.