The detection of high-energy neutrinos in the EeV range requires new detection techniques to cope with the small expected flux. The radio detection method, utilizing Askaryan emission, can be used to detect these neutrinos in polar ice. The propagation of the radio pulses has to be modeled carefully to reconstruct the energy, direction, and flavor of the neutrino from the detected radio flashes. Here, we study the effect of birefringence in ice, which splits up the radio pulse into two orthogonal polarization components with slightly different propagation speeds. This provides useful signatures to determine the neutrino energy and is potentially important to determine the neutrino direction to degree precision. We calculated the effect of birefringence from first principles where the only free parameter is the dielectric tensor as a function of position. Our code, for the first time, can propagate full RF waveforms, taking interference due to changing polarization eigenvectors during propagation into account. The model is available open-source through the NuRa-dioMC framework. We compare our results to in-situ calibration data from the ARA and ARIANNA experiments and find good agreement for the available time delay measurements, improving the predictions significantly compared to previous studies. Finally, the implications and opportunities for neutrino detection are discussed.