Tuning model parameters for complex climate codes is a challenging task due to the expense of a single simulation and a large number of uncertain input parameters. Bayesian calibration typically requires infeasibly many model evaluations on-the-fly. To accelerate Bayesian inference, we rely on polynomial chaos (PC) surrogates that approximate model input-output maps efficiently. Furthermore, the calibration procedure is enhanced to incorporate model structural errors, often the dominant component of predictive uncertainty. Namely, we develop a general framework for a probabilistic representation of the structural error inside the model, followed by a simultaneous calibration of physical inputs and parameters representing the structural error. The resulting embedded model-error strategy conserves physical constraints, allows meaningful predictions of a full set of output quantities of interest (QoIs), disambiguates model error from data noise, and leads to predictions with attributable uncertainties. The developed workflow is implemented in UQ Toolkit (www.sandia.gov/uqtoolkit). Surrogate-enabled sensitivity analysis and parameter inference are demonstrated for ELM FATES given observations, as well as for a simplified land model within the OSCM SciDAC project.