Gravitational waves emitted by neutron star black hole (NSBH) mergers encode key properties of black holes and neutron stars - such as their size, maximum mass and spins. However, the high mass ratio and the presence of matter makes generating long and accurate waveforms from these systems hard with numerical relativity, and not much is known about the systematic uncertainties due to waveform modeling. We simulate gravitational waves from NSBH mergers by hybridizing numerical relativity waveforms with a recent numerical relativity surrogate. These signals are analyzed using a range of available waveform families, and statistical and systematic errors are reported. We find that at a network signal-to-noise ratio (SNR) of 30, statistical uncertainties are usually larger than systematic offsets, while at an SNR of 70 the two become comparable. The individual black hole and neutron star masses, as well as the mass ratios, are typically measured very precisely, though not always accurately at high SNR. At a SNR of 30 the neutron star tidal deformability can only be bound from above, while for louder sources it can be measured and constrained away from zero. All neutron stars in our simulations are non-spinning, but in no case we can constrain the neutron star spin to be smaller than ~ 0.4 (90% credible interval). Waveform families tuned specifically for NSBH signals typically yield the most accurate characterization of the source parameters compared to those tuned for binary black holes or binary neutron stars.