  We develop new strategies to build numerical relativity surrogate models for
eccentric binary black hole systems, which are expected to play an increasingly
important role in current and future gravitational-wave detectors. We introduce
a new surrogate waveform model, \texttt{NRSur2dq1Ecc}, using 47 nonspinning,
equal-mass waveforms with eccentricities up to $0.2$ when measured at a
reference time of $5500M$ before merger. This is the first waveform model that
is directly trained on eccentric numerical relativity simulations and does not
require that the binary circularizes before merger. The model includes the
$(2,2)$, $(3,2)$, and $(4,4)$ spin-weighted spherical harmonic modes. We also
build a final black hole model, \texttt{NRSur2dq1EccRemnant}, which models the
mass, and spin of the remnant black hole. We show that our waveform model can
accurately predict numerical relativity waveforms with mismatches $\approx
10^{-3}$, while the remnant model can recover the final mass and dimensionless
spin with absolute errors smaller than $\approx 5 \times 10^{-4}M$ and $\approx
2 \times10^{-3}$ respectively. We demonstrate that the waveform model can also
recover subtle effects like mode-mixing in the ringdown signal without any
special ad-hoc modeling steps. Finally, we show that despite being trained only
on equal-mass binaries, \texttt{NRSur2dq1Ecc} can be reasonably extended up to
mass ratio $q\approx3$ with mismatches $\simeq 10^{-2}$ for eccentricities
smaller than $\sim 0.05$ as measured at a reference time of $2000M$ before
merger. The methods developed here should prove useful in the building of
future eccentric surrogate models over larger regions of the parameter space.
