How accurately can we predict the Moon's phases over periods of thousands of years?
The Jet Propulsion Laboratory produces mathematical models for the orbital motions of the Sun, Moon and planets. These are the models that are used to navigate spacecraft such as Voyager, Galileo and Cassini over billions of kilometres to rendezvous with the planets of the outer solar system.
These models are used, indirectly, as the basis for the Blue Moon calculator. They are very accurate indeed between the years 1800 and 2100, but like any mathematical models based on real data which contains measurement errors, their accuracy will decrease as the date is extended further into the past or future.
In particular, as the Moon slowly spirals away from the Earth, its angular speed is decreasing. This is called tidal deceleration. It amounts to about 25 arc-seconds per century per century, but we don't know the exact figure, and so we can't calculate the longitude of the Moon in its orbit with absolute certainty far into the future. An error of just 1 arc-second per century per century amounts to an uncertainty of more than an hour and a half in the time of Full Moon five thousand years ago or five thousand years hence.
There's another problem: the Earth itself, whose rate of rotation is not constant.
Models of the orbital motion of the Sun and Moon use a timescale called Dynamical Time, but civil timescales based on Greenwich Mean Time are defined in terms of the Earth's rotation. The Earth is slowing down, like an old-fashioned clock that needs to be wound up. As a result, GMT and Dynamical Time are drifting apart.
The Blue Moon calculator and the Moon Phase calculator both give dates and times expressed in GMT, using an approximation to the difference between GMT and Dynamical Time. This approximation is a quadratic expression for dates after the year 2000, and that will be reasonably accurate for a few decades into the 21st century.
By the year 7000, the quadratic formula predicts that the difference between the two timescales will be 24 hours, but this has no physical meaning and it should be taken with a pinch of salt.
The discussion of the frequency of Blue Moons uses Dynamical Time.