The anharmonic phonon behavior in zirconium hydrides and deuterides, including ϵ-ZrH2, γ-ZrH, and γ-ZrD, has been investigated from aspects of inelastic neutron scattering (INS) and lattice dynamics calculations within the framework of density functional theory (DFT). The harmonic model failed to reproduce the spectral features observed in the experimental data, indicating the existence of anharmonicity in those materials and the necessity of further explanations. Here, we present a detailed study on the anharmonicity in zirconium hydrides/deuterides by exploring the 2D potential energy surface of hydrogen/deuterium atoms and solving the corresponding 2D single-particle Schrödinger equation to obtain the eigenfrequencies, which are then convoluted with the instrument resolution. The convoluted INS spectra qualitatively describe the anharmonic peaks in the experimental INS spectra and demonstrate that the anharmonicity originates from the deviations of hydrogen potentials from quadratic behavior in certain directions; the effects are apparent for the higher-order excited vibrational states, but small for the ground and first excited states.