We propose two ensembles of geomagnetic field models spanning the last three millennia: COV-ARCH is calculated using all available archeological artefacts and volcanic lava flows; lake and marine sediment records are added to this data set to build COV-LAKE. Given the sparse distribution of archeomagnetic observations and their associated large uncertainties, the recovery of magnetic field changes from such data is an ill-posed inverse problem that requires assuming some prior information. This is usually performed by imposing arbitrary regularizations in space and time. Instead, we construct the prior knowledge entering the objective function to be minimized from spatial and temporal statistics of the geomagnetic field, as available from satellites, ground-based observatories and paleomagnetic measurements, and validated by numerical geodynamo simulations. Our approach relies on the projection of model coefficients onto temporal cross-covariance functions. We show with synthetic experiments that the dispersion within the ensemble of solutions provides a coherent measure of the model uncertainties. Gauss coefficients inverted from geophysical records compare satisfactorily with those deduced from the independent database built upon historical and observatory records. A posteriori model errors are reduced when incorporating sediment records; they logically increase towards decreasing length-scales, indicating that a partial information is available up to a spherical degree 4-5. Such models and their associated uncertainties are suited to be used as observations in geomagnetic data assimilation studies. Our results advocate for an approximately constant dipole decay since ≈1700 AD, preceded by an era (going back to 1000 AD) where the dipole trend is weak, possibly slightly positive. We observe in both hemispheres, at both low-and high-latitude, persistent patches over the past 3000 yr. We also confirm a westward drift of flux patches at the core-mantle boundary at a speed of about 0.20 to 0.25 • yr −1 . Despite the sparse data distribution in the southern hemisphere, the South Atlantic Anomaly appears in both ensembles of models around 1800 AD. A similar low-intensity event seems to have appeared below the Indian Ocean over 600-1400 AD. Both global models are in general good agreement with regional master curves, though filtering out some of the centennial oscillations.