[1] Inferring the core dynamics responsible for the observed geomagnetic secular variation requires knowledge of the magnetic field at the core-mantle boundary together with its associated model covariances. However, most currently available field models have been built using regularization conditions, which force the expansions in the spatial and time domains to converge but also hinder the calculation of reliable second-order statistics. To tackle this issue, we propose a stochastic approach that integrates, through time covariance functions, some prior information on the time evolution of the geomagnetic field. We consider the time series of spherical harmonic coefficients as realizations of a continuous and differentiable stochastic process. Our specific choice of process, such that it is not twice differentiable, mainly relies on two properties of magnetic observatory records (time spectra, existence of geomagnetic jerks). In addition, the required characteristic times for the low degree coefficients are obtained from available models of the magnetic field and its secular variation based on satellite data. We construct the new family COV-OBS of field models spanning the observatory and satellite era of 1840-2010. These models include the external dipole and permit sharper time changes of the internal field compared to previous regularized reconstructions. The a posteriori covariance matrix displays correlations in both space and time, which should be accounted for through the secular variation error model in core flow inversions and geomagnetic data assimilation studies.