We consider propagation of high energy earth-skimming taus produced in interactions of astrophysical tau neutrinos. For astrophysical tau neutrinos we take generic power-law flux, $E^{-2}$ and the cosmogenic flux initiated by the protons. We calculate tau energy loss in several approaches, such as dipole models and the phenomenological approach in which parameterization of the $F_2$ is used. We evaluate the tau neutrino charged-current cross section using the same approaches for consistency. We find that uncertainty in the neutrino cross section and in the tau energy loss partially compensate giving very small theoretical uncertainty in the emerging tau flux for distances ranging from $2$ km to $100$ km and for the energy range between $10^6$ GeV and $10^{11}$ GeV, focusing on energies above $10^8$ GeV. When we consider uncertainties in the neutrino cross section, inelasticity in neutrino interactions and the tau energy loss, which are not correlated, i.e. they are not all calculated in the same approach, theoretical uncertainty ranges from about $30\%$ and $60 \%$ at $10^8$ GeV to about factors of 3.3 and 3.8 at $10^{11}$ GeV for the $E^{-2}$ flux and the cosmogenic flux, respectively, for the distance of 10 km rock. The spread in predictions significantly increases for much larger distances, e.g., $\sim 1,000$ km. Most of the uncertainty comes from the treatment of photonuclear interactions of the tau in transit through large distances. We also consider Monte Carlo calculation of the tau propagation and we find that the result for the emerging tau flux is in agreement with the result obtained using analytic approach. Our results are relevant to several experiments that are looking for skimming astrophysical taus, such as the Pierre Auger Observatory, HAWC and Ashra. We evaluate the aperture for the Auger and discuss briefly application to the the other two experiments.