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 1 s 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 = 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 10 Hz. The surrogate has maximum errors of 3.8% in amplitude (0.04% excluding the last 100M before merger) and 0.043 radians in phase. This leads to typical mismatches of 10 −5 -10 −4 for Advanced LIGO depending on the component masses, with a worst case match of 7 × 10 −4 when both stars have masses ≥ 2M . The version implemented in the LIGO Algorithm Library takes ∼ 0.07 s to evaluate for a starting frequency of 30 Hz and ∼ 0.8 s for a starting frequency of 10 Hz, resulting in a speed-up factor of O(10 3 ) 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.