Paleoclimate reconstructions aim to recreate past climates and are critical in assessing how modern day temperatures (and other climate variables) are anomalous in a millennial context. Most of the methods available in the literature give consolidated estimates of climate variables over the past millennium and beyond. It is widely recognized that changes in climate differ by spatial location; for example, estimated changes in temperature, and variability thereof, between the tropics and the polar regions are different. Obtaining a detailed understanding of such phenomena requires developing methodology for paleoclimate reconstructions that are spatially disaggregated, with uncertainty quantifications of estimated past climate. Spatially disaggregate paleoclimate reconstructions are fundamentally high dimensional statistical problems, with unique challenges imposed by the geosciences context. In particular, paleoclimate reconstruction rely on relating climate variables to proxy variables (such as tree rings, ice cores etc..). The project will lay the theoretical foundation for high dimensional paleoclimate reconstructions using modern day statistical methods. As a concrete application the project will reconstruct global past climates at a spatially disaggregated levels for the past millennium, and also attach confidence statements to these reconstructions. The proposed work entails a collaboration between the mathematical and geo-sciences to solve a scientific problem at the interface of both fields. The methods developed during the project also have broad applications in other fields, such as genomics, where relationships between genes in a high dimensional context are often studied. Climate field reconstructions are inherently multivariate inference problems (e.g., spatial data on a large grid are required), rely on noisy input data, and are often so high-dimensional that the data dimension is close to or exceeds the sample size, resulting in ill-posed or rank-deficient estimation problems. In this context high dimensional mean and covariance estimation is often at the center of the inferential problem. Furthermore, meaningful solutions to such problems require a reliable knowledge of the uncertainty in estimated model parameters. The need for a rigorous quantification of uncertainties has recently spurred much interest in Bayesian methods for climate reconstruction, usually based on Markov Chain Monte Carlo (MCMC) sampling techniques. The latter pose two major issues, as: (1) there is no guarantee that posterior samples are generated from the required distribution (convergence issue); (2) they rely on computationally-heavy algorithms which limit their applicability (efficiency, scalability and applicability). We propose to overcome both limitations by using a flexible but high dimensional Bayesian approach that leads to closed-form solutions for posterior quantities, hence alleviating much of the computational burden. Convergence issues will be addressed first, and dimensionality reduction will subsequently be implemented by exploiting the rich theory of Markov Random fields. The outcome of the theoretical component of the project will be the construction of novel and sound, efficient high-dimensional Bayesian algorithms tailor-made for climate field reconstruction problems. These new statistical tools will then be applied to the reconstruction of global and regional temperature fields from heterogeneous geological proxies (tree rings, ice cores, speleothems, corals, sediments) over the Common Era, and new analysis of instrumental surface temperature and sea-level pressure. The solutions will be accompanied by credible intervals, and promises to yield new insights into natural climate variability at global scales.