Archaeomagnetic observations are known for being most challenging when it comes to inferring the geomagnetic field at their origin, due mainly to their poor global coverage and substantial age uncertainties. To supplement the information provided by the sparse archaeomagnetic dataset, we consider the extra information on the magnetic field given by numerical simulations of the geodynamo. Ideally, both data and numerical models would be connected optimally through data assimilation methods (for instance using a Kalman Filter), in order to estimate the magnetic field, fluid flow and buoyancy in the Earth's core. However, numerical models are substantially nonlinear and the archaeomagnetic data nonlinearly connected to the magnetic field, spoiling the optimality of the solution. We thus resort to an Ensemble Kalman Filter in order to deal partly with the nonlinearity of the model, and we apply a linearization of the observation operator in order to assimilate the data into the numerical model. As a first test, we apply the assimilation method to archaeomagnetic data, producing estimates of the magnetic field at the core surface through the last three millennia. Further, to access the sensitivity of the scheme to the prior information, we make use of different dynamo simulations, of varying complexity, as sources of prior information. Upon validation, we aim at applying the data assimilation algorithm sequentially. We expect that each analysis of the magnetic field given the available data will lead to a propagation of the new information throughout the core by the numerical model, making it possible to better understand the state of the core through the last millennia.