Global field reconstructions of the past are a key tool for understanding the dynamics of the Earth's magnetic field and the underlying processes in the Earth's core (e.g., Constable & Korte, 2015). This includes studying the evolution of field features, such as dipole decay, the South Atlantic Anomaly (SAA) and flux patches (Hartmann & Pacca, 2009;Jackson & Finlay, 2015). In the past, several techniques for constructing global field models have been developed and employed. Truncated spherical harmonics (SH) in the spatial domain combined with spline interpolation in time are widely used (Jackson et al., 2000;Korte et al., 2009). In the eighties, Constable and Parker (1988) first proposed using Gaussian processes to model the field dynamics, but until recently, the technique had not been applied to global field modeling. Only in the last years, statistical methods implementing this approach have been suggested (Hellio & Gillet, 2018;Nilsson & Suttie, 2021).Even though Bloxham and Jackson (1992) already discussed the accurate assessment of uncertainties, most following studies did not pursue the suggested analytical approaches, and either did not report uncertainties at