Climate and land surface model intercomparison projects give a qualitative indication of the range of behaviors one may expect from a group of 'realistic' models. As stated within the IPCC report Climate Change 2001: The Scientific Basis, there needs to be a 'more comprehensive and systematic system to of model analysis and diagnosis, and a Monte Carlo approach to model uncertainties associated with parameterizations and initial conditions'. One of the reasons these sorts of calculations have been limited is the computational cost of evaluating the uncertainty associated with any given choice of model parameters or the uncertainty associated with any given combination of model parameters. We present a new method based on Bayesian Stochastic Inversion that can efficiently determine the multi-dimensional probability density functions and error covariances. When optimized, BSI is several orders of magnitude more efficient than analyses based on the Monte Carlo approach. We present the results of a study using a surrogate climate system model to estimate how noise, degree of non-linearity, and number of parameters affect the cost of the analysis. We also illustrate the power of this approach in determining the optimal parameters and estimates of uncertainty within off-line predictions of land surface energy fluxes at the Cabauw data-site, Netherlands, using the CHASM (CHAmeleon Surface Model) land surface model.