  Numerical relativity simulations provide the most precise templates for the
gravitational waves produced by binary black hole mergers. However, many of
these simulations use an incomplete waveform extraction technique --
extrapolation -- that fails to capture important physics, such as gravitational
memory effects. Cauchy-characteristic evolution (CCE), by contrast, is a much
more physically accurate extraction procedure that fully evolves Einstein's
equations to future null infinity and accurately captures the expected physics.
In this work, we present a new surrogate model, NRHybSur3dq8$\_$CCE, built from
CCE waveforms that have been mapped to the post-Newtonian (PN) BMS frame and
then hybridized with PN and effective one-body (EOB) waveforms. This model is
trained on 102 waveforms with mass ratios $q\leq8$ and aligned spins
$\chi_{1z}, \, \chi_{2z} \in \left[-0.8, 0.8\right]$. The model spans the
entire LIGO-Virgo-KAGRA (LVK) frequency band (with
$f_{\text{low}}=20\text{Hz}$) for total masses $M\gtrsim2.25M_{\odot}$ and
includes the $\ell\leq4$ and $(\ell,m)=(5,5)$ spin-weight $-2$ spherical
harmonic modes, but not the $(3,1)$, $(4,2)$ or $(4,1)$ modes. We find that
NRHybSur3dq8$\_$CCE can accurately reproduce the training waveforms with
mismatches $\lesssim2\times10^{-4}$ for total masses $2.25M_{\odot}\leq
M\leq300M_{\odot}$ and can, for a modest degree of extrapolation, capably model
outside of its training region. Most importantly, unlike previous waveform
models, the new surrogate model successfully captures memory effects.
