We address the problem of reconstructing a global field from proxy data with sparse spatial sampling such as the MARGO (multiproxy approach for the reconstruction of the glacial ocean surface) SST (sea-surface temperature) and d18Oc (oxygen-18/oxygen-16 isotope ratio preserved in fossil carbonate shells of planktic foraminifera) data. To this end, we propose to assimilate these data into coupled climate models by adjusting some of their parameters and optimizing the fit. In particular, we suggest to combine a forward model and an objective function that quantifies the misfit to the data. Because of their computational efficiency, earth system models of intermediate complexity are particularly well-suited for this purpose. We used one such model (the University of Victoria EarthSystem Climate Model) and carried out a series of sensitivity experiments by varying a single model parameter through changing the atmospheric CO2 concentration. The unanalyzed World Ocean Atlas SST and the observed sea-ice concentration served as present-day targets. The sparse data coverage as implied by the locations of 756 ocean sediment cores from the MARGO SST database was indeed suf cient to determine the best fit. As anticipated, it turned out to be the 365 ppm experiment. We also found that the 200 ppm experiment came surprisingly close to what is commonly expected for the Last Glacial Maximum ocean circulation. Our strategy has a number of advantages over more traditional mapping methods, e.g., there is no need to force the results of different proxies into a single map, because they can be compared to the model output one at a time, properly taking into account the different seasons of plankton growth or varying depth habitats. It can be extended to more model parameters and even be automated.