Predicting the properties of the matter ejected during and after a neutron star merger is crucial to our ability to use electromagnetic observations of these mergers to constrain the masses of the neutron stars, the equation of state of dense matter, and the role of neutron star mergers in the enrichment of the Universe in heavy elements. Our ability to reliably provide such predictions is however limited by a broad range of factors, including the finite resolution of numerical simulations, their treatment of magnetic fields, neutrinos, and neutrino-matter interactions, and the approximate modeling of the equation of state of dense matter. To disentangle some of these uncertainties, in this manuscript we focus narrowly on the impact that two small variations in the setup of numerical simulations can have on the matter ejected during the merger. Specifically, we perform 4 simulations: one at higher residual eccentricity (e ∼ 0.01), two with reduced eccentricity (e ∼ 0.003–0.004), and slightly different equation of state implementations, and a fourth at higher resolution to estimate numerical errors. We perform these simulations for a single physical configuration, a 1.3M⊙−1.4M⊙ binary neutron star system using the SFHo equation of state, ignoring magnetic fields and neutrino effects. We find that a residual eccentricity e ∼ 0.01, as measured ∼4–6 orbits before merger, causes O(25–30%) relative changes in the amount of ejected mass. These differences in mass ejection appear mainly due to changes in the amount of matter ejected as a result of core bounces during merger. We note that O(1%) residual eccentricities have regularly been used in binary neutron star merger simulations as proxy for circular binaries, potentially creating an additional source of error in predictions for the mass of the dynamical ejecta. On the other hand, we find that simulations performed with slightly different implementations of the same equation of state are consistent within estimated numerical errors.