  Gravitational waves from compact binary coalescence sources can be decomposed
into spherical-harmonic multipoles, the dominant being the quadrupole ($\ell=2,
m=\pm2$) modes. The contribution of sub-dominant modes towards total signal
power increases with increasing binary mass ratio and source inclination to the
detector. It is well-known that in these cases neglecting higher modes could
lead to measurement biases, but these have not yet been quantified with a
higher-mode model that includes spin effects. In this study, we use the
multi-mode aligned-spin phenomenological waveform model IMRPhenomHM to
investigate the effects of including multi-mode content in estimating source
parameters and contrast the results with using a quadrupole-only model
(IMRPhenomD). We use as sources IMRPhenomHM and hybrid EOB-NR waveforms over a
range of mass-ratio and inclination combinations, and recover the parameters
with IMRPhenomHM and IMRPhenomD. These allow us to quantify the accuracy of
parameter measurements using a multi-mode model, the biases incurred when using
a quadrupole-only model to recover full (multi-mode) signals, and the
systematic errors in the IMRPhenomHM model. We see that the parameters
recovered by multi-mode templates are more precise for all non-zero
inclinations as compared to quadrupole templates. For multi-mode injections,
IMRPhenomD recovers biased parameters for non-zero inclinations with lower
likelihood while IMRPhenomHM recovered parameters are accurate for most cases,
and if a bias exists, it can be explained as a combined effect of observational
priors and (in the case of hybrid-NR signals) waveform inaccuracies. For cases
where IMRPhenomHM recovers biased parameters, the bias is always smaller than
the corresponding IMRPhenomD recovery, and we conclude that IMRPhenomHM will be
sufficiently accurate to allow unbiased measurements for most GW observations.
