Most natural interpolation to simulate increasing distance?
One typically interpolates going up, 16 $\times$16 to 32 $\times$ 32.  Going down is low pass filtering (anti aliasing) and decimation.  As you have acknowledged, this isn’t the same as the effect of distance.  As long you understand that your original image was captured with a camera with a specific MTF, and sensor resolution, and are going to neglect scattering (which would require something like MODTRAN),  your original image has been acquired with some angular resolution, $\theta_0$, your image has a linear spatial resolution $\Delta_0$ such that $$
\Delta_0 = r_0 \theta_0
$$ where $r_0$ is the distance that the image was acquired at.  At some greater distance $r_g$, $\theta_0$ can be considered the same so,
$$ \Delta_g= \theta_0 r_g. $$
Given some coordinate system associated with the original system, you could use something like a Gaussian blur with a standard deviation related to $\Delta_g$, and resample. This doesn’t really help with a large depth of field.