The interpretation of cosmological observations relies on a notion of average Universe, which is usually taken as the homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) model. However, inhomogeneities may statistically bias the observational averages with respect to FLRW, notably for distance measurements, due to a number of effects such as gravitational lensing and redshift perturbations. In this article, we review the main known theoretical results on average distance measures in cosmology, based on second-order perturbation theory, and fill some of their gaps. We then comprehensively test these theoretical predictions against ray tracing in a high-resolution dark-matter N-body simulation. This method allows us to describe the effect of small-scale inhomogeneities deep into the non-linear regime of structure formation, on light propagation up to z = 10. We find that numerical results are in remarkably good agreement with theoretical predictions, in the limit of super-sample variance. No unexpectedly large bias originates from the very small scales, whose effect is fully encoded in the non-linear power spectrum. Specifically, the directional average of the inverse amplification and the source-averaged amplification are compatible with unity; the change in area of surfaces of constant cosmic time is compatible with zero; the biases on other distance measures, which can reach ∼ 10 −3 at high redshift, are well understood. As a side product, we also confront the predictions of the recent finite-beam formalism with numerical data, and find excellent agreement.