As current gravitational wave (GW) detectors increase in sensitivity, and particularly as new instruments are being planned, there is the possibility that ground-based GW detectors will observe GWs from highly eccentric neutron star binaries. We present the first detailed study of highly eccentric BNS systems with full (3+1)D numerical relativity simulations using consistent initial conditions, i.e., setups which are in agreement with the Einstein equations and with the equations of general relativistic hydrodynamics in equilibrium. Overall, our simulations cover two different equations of state (EOSs), two different spin configurations, and three to four different initial eccentricities for each pairing of EOS and spin. We extract from the simulated waveforms the frequency of the f-mode oscillations induced during close encounters before the merger of the two stars. The extracted frequency is in good agreement with f-mode oscillations of individual stars for the irrotational cases, which allows an independent measure of the supranuclear equation of state not accessible for binaries on quasicircular orbits. The energy stored in these f-mode oscillations can be as large as 10−3 M⊙∼1051 erg, even with a soft EOS. In order to estimate the stored energy, we also examine the effects of mode mixing due to the stars’ offset from the origin on the f-mode contribution to the GW signal. While in general (eccentric) neutron star mergers produce bright electromagnetic counterparts, we find that for the considered cases with fixed initial separation the luminosity decreases when the eccentricity becomes too large, due to a decrease of the ejecta mass. Finally, the use of consistent initial configurations also allows us to produce high-quality waveforms for different eccentricities which can be used as a test bed for waveform model development of highly eccentric binary neutron star systems.S. V. C. was supported by the DFG Research Training Group 1523/2 “Quantum and Gravitational Fields.” T. D. acknowledges support by the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 749145, BNSmergers. N. K. J.-M. acknowledges support from STFC Consolidator Grant No. ST/L000636/1. B. B. was supported by DFG Grant No. BR 2176/5-1. W. T. was supported by the National Science Foundation under Grant No. PHY-1707227. Also, this work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 690904. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation. Computations were performed on the supercomputer SuperMUC at the LRZ (Munich) under the project number pr48pu and on the ARA cluster of the University of Jena