Over the past year, a handful of new gravitational wave models have been developed to include multiple harmonic modes thereby enabling for the first time fully Bayesian inference studies including higher modes to be performed. Using one recently-developed numerical relativity surrogate model, NRHybSur3dq8, we investigate the importance of higher modes on parameter inference of coalescing massive binary black holes. We focus on examples relevant to the current three-detector network of observatories, with detector-frame masses and signal amplitudes consistent with plausible candidates for the next few observing runs. We show that for such systems the higher mode content will be important for interpreting coalescing binary black holes, reducing systematic bias, and computing properties of the remnant object. For asymmetric binaries with mass ratios q>1, higher modes are critical, and their omission usually produces substantial parameter biases. Even for comparable-mass binaries and at low signal amplitude, the omission of higher modes can influence posterior probability distributions. We discuss the impact of our results on source population inference and self-consistency tests of general relativity.