  Gravitational-wave detectors have begun to observe coalescences of heavy
black holes at a consistent pace for the past few years. Accurate models of
gravitational waveforms are essential for unbiased and precise estimation of
source parameters. Recently developed surrogate models based on high-accuracy
numerical relativity (NR) simulations are ideal for constraining physical
parameters of heavy black hole merger events. In this paper, we first
demonstrate the viability of these multi-modal surrogates as reliable parameter
estimation tools. We show that NR surrogates can extract additional information
from GW data that is inaccessible to traditional models, by analyzing a set of
synthetic signals with the NR surrogate and other approximants.
  We also consider the case of two of the earliest binary black holes detected
by the LIGO observatories: GW150914 and GW170104. We reanalyze their data with
fully-precessing NR-surrogate templates and freely provide the resulting
posterior samples as supplemental material. We find that our refined analysis
is able to extract information from sub-dominant GW harmonics in data, and
therefore better resolve the degeneracy in measuring source luminosity distance
and orbital inclination for both events. We estimate the sources of both events
to be 20-25% further away than was previously estimated. Our analyses also
constrain their orbital orientation more tightly around face-on or face-off
configurations than before. Additionally, for GW150914 we constrain the
effective inspiral spin more tightly around zero. This work is one of the first
to unambiguously extract sub-dominant GW mode information from real events. It
is also a first step toward eliminating the approximations used in
semi-analytic waveform models from GW parameter estimation. It also motivates
that NR surrogates be extended to cover more of the binary black hole parameter
space.
