  Simulating a binary black hole (BBH) coalescence by solving Einstein's
equations is computationally expensive, requiring days to months of
supercomputing time. Using reduced order modeling techniques, we construct an
accurate surrogate model, which is evaluated in a millisecond to a second, for
numerical relativity (NR) waveforms from non-spinning BBH coalescences with
mass ratios in $[1, 10]$ and durations corresponding to about $15$ orbits
before merger. We assess the model's uncertainty and show that our modeling
strategy predicts NR waveforms {\em not} used for the surrogate's training with
errors nearly as small as the numerical error of the NR code. Our model
includes all spherical-harmonic ${}_{-2}Y_{\ell m}$ waveform modes resolved by
the NR code up to $\ell=8.$ We compare our surrogate model to Effective One
Body waveforms from $50$-$300 M_\odot$ for advanced LIGO detectors and find
that the surrogate is always more faithful (by at least an order of magnitude
in most cases).
