  Coalescing binary black holes are among the primary science targets for
second generation ground-based gravitational wave (GW) detectors. Reliable GW
models are central to detection of such systems and subsequent parameter
estimation. This paper performs a comprehensive analysis of the accuracy of
recent waveform models for binary black holes with aligned spins, utilizing a
new set of $84$ high-accuracy numerical relativity simulations. Our analysis
covers comparable mass binaries ($1\le m_1/m_2\le 3$), and samples
independently both black hole spins up to dimensionless spin-magnitude of $0.9$
for equal-mass binaries and $0.85$ for unequal mass binaries. Furthermore, we
focus on the high-mass regime (total mass $\gtrsim 50M_\odot$). The two most
recent waveform models considered (PhenomD and SEOBNRv2) both perform very well
for signal detection, losing less than 0.5\% of the recoverable signal-to-noise
ratio $\rho$, except that SEOBNRv2's efficiency drops mildly for both black
hole spins aligned with large magnitude. For parameter estimation, modeling
inaccuracies of SEOBNRv2 are found to be smaller than systematic uncertainties
for moderately strong GW events up to roughly $\rho\lesssim 15$. PhenomD's
modeling errors are found to be smaller than SEOBNRv2's, and are generally
irrelevant for $\rho\lesssim 20$. Both models' accuracy deteriorates with
increased mass-ratio, and when at least one black hole spin is large and
aligned. The SEOBNRv2 model shows a pronounced disagreement with the numerical
relativity simulation in the merger phase, for unequal masses and
simultaneously both black hole spins very large and aligned. Two older waveform
models (PhenomC and SEOBNRv1) are found to be distinctly less accurate than the
more recent PhenomD and SEOBNRv2 models. Finally, we quantify the bias expected
from all GW models during parameter estimation for recovery of binary's masses
and spins.
