For the first time, we construct an inspiral-merger-ringdown waveform model within the effective-one-body formalism for spinning, nonprecessing binary black holes that includes gravitational modes beyond the dominant (?,|m|)=(2,2) mode, specifically (?,|m|)=(2,1),(3,3),(4,4),(5,5). Our multipolar waveform model incorporates recent (resummed) post-Newtonian results for the inspiral and information from 157 numerical-relativity simulations, and 13 waveforms from black-hole perturbation theory for the (plunge-)merger and ringdown. We quantify the improvement in accuracy when including higher-order modes by computing the faithfulness of the waveform model against the numerical-relativity waveforms used to construct the model. We define the faithfulness as the match maximized over time, phase of arrival, gravitational-wave polarization and sky position of the waveform model, and averaged over binary orientation, gravitational-wave polarization and sky position of the numerical-relativity waveform. When the waveform model contains only the (2,2) mode, we find that the averaged faithfulness to numerical-relativity waveforms containing all modes with ??5 ranges from 90% to 99.9% for binaries with total mass 20–200 M? (using the Advanced LIGO’s design noise curve). By contrast, when the (2,1), (3,3), (4,4), (5,5) modes are also included in the model, the faithfulness improves to 99% for all but four configurations in the numerical-relativity catalog, for which the faithfulness is greater than 98.5%. Starting from the complete inspiral-merger-ringdown model, we develop also a (stand-alone) waveform model for the merger-ringdown signal, calibrated to numerical-relativity waveforms, which can be used to measure multiple quasi-normal modes. The multipolar waveform model can be extended to include spin-precessional effects, and will be employed in upcoming observing runs of Advanced LIGO and Virgo.