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.