Numerical geodynamo simulations capture several features of the spatial and temporal geomagnetic field variability on historical and Holocene timescales. However, a recent analysis questioned the ability of these numerical models to comply with long-term palaeomagnetic field behaviour. Analysing a suite of 50 geodynamo models, we present here the first numerical simulations known to reproduce the salient aspects of the palaeosecular variation and time-averaged field behaviour since 10 Ma. We find that the simulated field characteristics covary with the relative dipole field strength at the core-mantle boundary (dipolarity). Only models dominantly driven by compositional convection, with an Ekman number (ratio of viscous to Coriolis forces) lower than $10^{-3}$ and a dipolarity in the range $0.34-0.56$ can capture the observed palaeomagnetic field behaviour. This dipolarity range agrees well with state-of-the-art statistical field models and represent a testable prediction for next generation global palaeomagnetic field model reconstructions.