The focus of this thesis is to develop a general methodology to obtain high-resolution spatial-temporal forecasts of Sea Surface Temperature (SST) using an ensembles of general circulation model (GCM) output and historical records as the major driving force. As a case study, we consider Sea Surface Temperature (SST) in the North Pacific Ocean. We use two ensembles of different GCM simulation output, made available in the 4th Assessment Report of the Intergovernmental Panel on Climate Change: one corresponds to 20th century forcing conditions and the other corresponds to the A1B emissions scenario for the 21st century. Given a representation of the SST spatio-temporal fields based on a common set of empirical orthogonal functions (EOFs), we use a hierarchical Bayesian model for the EOF coefficients to estimate a baseline and a set of model discrepancies. These components are all time-varying. The model enables us to extract relevant temporal patterns of variability from both the observations and simulations and obtain common patterns from all eighteen series. This is used to obtain unified 21st century forecasts of relevant oceanic indexes as well as whole fields of forecast North Pacific SST. The unified forecast captures large longterm oceanic behavior, however the coarse resolution prevents us from capturing coastal behaviors. We use the unified forecast to model high resolution SST by establishing a link between large and small scale variability using statistical downscaling techniques. Using a combination of a discrete process convolution and a dynamic linear model, we obtain a smooth high-resolution forecast of SST fields off the coast of California. To model the high resolution data faster and efficiently, we developed and implement a parallel version of the forward filtering backwards sampling algorithm. We finish the work with remarks on the model results and address future avenues this work can take.