Model structural error is often the dominant component of predictive uncertainty budget in climate models. 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 approach is further enhanced to include a spatio-temporal model surrogate with Karhunen-Loeve and polynomial chaos representations, providing dimensionality reduction with quantifiable uncertainty, and augmenting the predictive variance. The developed workflow is implemented in UQ Toolkit (www.sandia.gov/uqtoolkit). The method is demonstrated for E3SM (Energy Exascale Earth System Model) land model calibration given FLUXNET observations, highlighting the need for burdening physical parameters with stochasticity due to forcing factors.