The use of Earth System Models in many-query applications, such as uncertainty quantification and optimization, is limited by their high computational cost. Surrogates that approximate the input-output maps are typically constructed and employed for computationally intensive tasks. We will demonstrate methods to incorporate additional stochasticity into these surrogates to characterize uncertainties arising from these models' intrinsic stochasticity. Such generative, parametric-stochastic surrogates will then be employed in a model calibration loop in which both physical parameters and parameters associated with structural model errors are calibrated to fit observational data. With such embedded model-error statistical representations, the calibrated model predictions remain consistent with the Earth System Model’s physical constraints, thanks to the non-intrusive nature of the model-error embedding. Besides, the calibrated model is endowed with predictive uncertainty, which can be decomposed into components such as intrinsic stochasticity, observational data noise, parametric uncertainty, surrogate errors, and model structural errors. The methods will be demonstrated on various forms of the Energy Exascale Earth System Model (E3SM), including the land component and the analysis of the quasi-biennial oscillation.