@PM2Ring's answer does a great job of explaining that folks who do careful numerical integration keep a large number of digits beyond the final accuracy because if they didn't, after tens of thousands or even millions of integration steps computer round-off error would accumulate and threaten to dominate as a source of error.
Since I don't yet know how many steps they use per year of simulated time, I've just asked
Further, for most computer architectures and libraries we don't have the option to use exactly the number of digits of accuracy we need. Instead there are standards that define 32, 64, and 128 bit binary representations of floating numbers for example, and computers and their numerical processors are built around these.
It's dangerous to write anything about this topic without an expert dropping by and correcting this or that because numerical precision is such a critical issue, but to illustrate I'll type
numpy.finfo(float) into Python and get for example
finfo(resolution=1e-15, min=-1.7976931348623157e+308, max=1.7976931348623157e+308, dtype=float64)
We can see that the computer uses 64 bits for mantissa, exponent, and two signs. Sixteen decimal places (your number) eats up roughly 53 bits out of 64, the rest go to the exponent and two signs.
It's pretty common to share raw output of the numerical integration by displaying all the digits. Since solar system distances when expressed in kilometers can be shown without needing scientific notation, you end up with sixteen digits displayed.
The underlying question
How on Earth (no pun intended) do astronomers make measurements of planetary position accurate to... metres?
As pointed out, timing pulses of light reflected from the Moon doesn't directly measure the center of the moon, but measuring round trip delays and doppler shifts of microwaves from Earth to spacecraft orbiting solar system bodies for weeks, months, or even years does allow one to eventually get the distances to the gravitational centers of these bodies to meters.
That's done by simulating, varying the inputted masses and finding the best match to the data for the microwave delays and Doppler shifts collected.
We can read more about the ephemerides (pronounciation) used in JPL Horizons On-Line Ephemeris System
Since August 2013, Horizons has been using ephemeris DE431. During the week of 12 April 2021, the Horizons ephemeris system was updated to replace the DE430/431 planetary ephemeris, used since 2013, with the new DE440/441 solution. The new DE440/441 general-purpose planetary solution includes seven additional years of ground and space-based astrometric data, data calibrations, and dynamical model improvements, most significantly involving Jupiter, Saturn, Pluto, and the Kuiper Belt. Inclusion of 30 new Kuiper-belt masses, and the Kuiper Belt ring mass, results in a time-varying shift of ~100 km in DE441's barycenter relative to DE431.
Let's look at The JPL Planetary and Lunar Ephemerides DE440 and DE441 You can read about the numerical integration and accuracy, but here I'll show how the final results match to the measure data for orbiters around solar system bodies.
You can see that the residual errors are indeed of the order of a few meters. This is a result of simulating the many gravitational forces between a large number of solar system bodies (planets, asteroids, moons, etc.) and some non-gravitational effects using equations similar to $F=ma$ except corrected to some order for relativistic effects.
For more on that see answers to
Not shown are any of the angular determinations for orbiting spacecraft made by VLBI.
Figure 4. Residuals of LLR ranges against DE440. The rms residual of the LLR ranges is about 20 cm for the early data and is about 1.3 cm for the recent data.
Figure 5. Residuals of the MESSENGER range data against DE440. The rms residual of the MESSENGER ranges is about 0.7 m.
Figure 6. Residuals of the Venus Express range data against DE440. The rms residual of the Venus Express ranges is about 8 m.
Figure 7. Residuals of the Mars orbiter range data against DE440. The rms residual of the MEX ranges is about 2 m and the rms residuals of the MGS, ODY, and MRO ranges are about 0.7 m.
Figure 9. Residuals of the Juno range data against DE440. The rms residual of the Juno ranges is about 13 m.