  We present a comprehensive comparison between numerical relativity (NR)
angular momentum fluxes at infinity and the corresponding quantity entering the
radiation reaction in TEOBResumS, an Effective-One-Body (EOB) waveform model
for nonprecessing coalescing black hole binaries on quasi-circular orbits. This
comparison prompted us to implement two changes in the model: (i) including
Next-to-Quasi-Circular corrections in the $\ell=m$, $\ell\leq 5$ multipoles
entering the radiation reaction and (ii) consequently updating the NR-informed
spin-orbital sector of the model. This yields a new waveform model that
presents a higher self-consistency between waveform and dynamics and an
improved agreement with NR simulations. We test the model computing the EOB/NR
unfaithfulness $\bar{F}_{\rm EOB/NR}$ over all 534 spin-aligned configurations
available through the Simulating eXtreme Spacetime catalog, notably using the
noise spectral density of Advanced LIGO, Einstein Telescope and Cosmic
Explorer, for total mass up to $500M_\odot$. We find that the maximum
unfaithfulness $\bar{F}^{\rm max}_{\rm EOB/NR}$ is mostly between $10^{-4}$ and
$10^{-3}$, and the performance progressively worsens up to $\sim 5\times
10^{-3}$ as the effective spin of the system is increased. We perform similar
analyses on the \SEOB{} model, that delivers $\bar{F}^{\rm max}_{\rm EOB/NR}$
values uniformly distributed versus effective spin and mostly between $10^{-3}$
and $10^{-2}$. We conclude that the improved TEOBResumS model already
represents a reliable and robust first step towards the development of highly
accurate waveform templates for third generation detectors.
