There are several Global Climate Models (GCM) reported by various countries to the Intergovernmental Panel on Climate Change (IPCC). Due to the varied nature of the GCM model assumptions, the future projections of the GCMs show high variability which makes it difficult to come up with confident projections into the future. Climate scientists combine these multiple GCM models to minimize the variability and the prediction error. Most of these model combinations are specifically for a location, or at a global scale. They do not consider regional or local smoothing (including the IPCC model). In this paper, we address this problem of combining multiple GCM model outputs with spatial smoothing as an important desired criterion. The problem formulation takes the form of multiple least squares regression for each geographic location with graph Laplacian based smoothing amongst the neighboring locations. Unlike the existing Laplacian regression frameworks, our formulation has both inner and outer products of the coefficient matrix, and has Sylvester equations as its special case. We discuss a few approaches to solve the problem, including a closed-form by solving a large linear system, as well as gradient descent methods which turn out to be more efficient. We establish the superiority of our approach in terms of model accuracy and smoothing compared to several popular baselines on real GCM climate datasets.