  Numerical relativity (NR) simulations provide the most accurate binary black
hole gravitational waveforms, but are prohibitively expensive for applications
such as parameter estimation. Surrogate models of NR waveforms have been shown
to be both fast and accurate. However, NR-based surrogate models are limited by
the training waveforms' length, which is typically about 20 orbits before
merger. We remedy this by hybridizing the NR waveforms using both
post-Newtonian and effective one body waveforms for the early inspiral. We
present NRHybSur3dq8, a surrogate model for hybridized nonprecessing numerical
relativity waveforms, that is valid for the entire LIGO band (starting at
$20~\text{Hz}$) for stellar mass binaries with total masses as low as
$2.25\,M_{\odot}$. We include the $\ell \leq 4$ and $(5,5)$ spin-weighted
spherical harmonic modes but not the $(4,1)$ or $(4,0)$ modes. This model has
been trained against hybridized waveforms based on 104 NR waveforms with mass
ratios $q\leq8$, and $|\chi_{1z}|,|\chi_{2z}| \leq 0.8$, where $\chi_{1z}$
($\chi_{2z}$) is the spin of the heavier (lighter) BH in the direction of
orbital angular momentum. The surrogate reproduces the hybrid waveforms
accurately, with mismatches $\lesssim 3\times10^{-4}$ over the mass range
$2.25M_{\odot} \leq M \leq 300 M_{\odot}$. At high masses
($M\gtrsim40M_{\odot}$), where the merger and ringdown are more prominent, we
show roughly two orders of magnitude improvement over existing waveform models.
We also show that the surrogate works well even when extrapolated outside its
training parameter space range, including at spins as large as 0.998. Finally,
we show that this model accurately reproduces the spheroidal-spherical mode
mixing present in the NR ringdown signal.
