  The interpretation of cosmological observations relies on a notion of average
Universe, which is usually taken as the homogeneous and isotropic
Friedmann-Lema\^itre-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 $\sim
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.
