Gravitational-wave observations of binary neutron star systems can provide information about the masses, spins, and structure of neutron stars. However, this requires accurate and computationally efficient waveform models that take <1s to evaluate for use in Bayesian parameter estimation codes that perform 10^7 - 10^8 waveform evaluations. We present a surrogate model of a nonspinning effective-one-body waveform model with l = 2, 3, and 4 tidal multipole moments that reproduces waveforms of binary neutron star numerical simulations up to merger. The surrogate is built from compact sets of effective-one-body waveform amplitude and phase data that each form a reduced basis. We find that 12 amplitude and 7 phase basis elements are sufficient to reconstruct any binary neutron star waveform with a starting frequency of 10Hz. The surrogate has maximum errors of 3.8% in amplitude (0.04% excluding the last 100M before merger) and 0.043 radians in phase. The version implemented in the LIGO Algorithm Library takes ~0.07s to evaluate for a starting frequency of 30Hz and ~0.8s for a starting frequency of 10Hz, resulting in a speed-up factor of ~10^3 - 10^4 relative to the original Matlab code. This allows parameter estimation codes to run in days to weeks rather than years, and we demonstrate this with a Nested Sampling run that recovers the masses and tidal parameters of a simulated binary neutron star system.