Computational models usually carry misspecification due to different parameterizations and assumptions. Ignoring such model errors can lead to poor predictive capability, even when high-quality data are used. It is thus crucial to quantify and propagate uncertainty due to model error, and to differentiate it from parametric uncertainty and data noise. Traditional approaches accommodate model error through discrepancy terms that are only available for calibration quantities, and generally do not preserve physical model constraints in subsequent predictions. The ability to extrapolate to other predictive quantities and to retain certain physical properties are often desirable and even required in many physical science and engineering applications. We develop a stochastically embedded model correction approach that enables these qualities, and perform Bayesian inference of the correction terms together with model parameters. Employing polynomial chaos expansions to represent the correction terms, the approach allows efficient quantification, propagation, and decomposition of uncertainty that includes contributions from data noise, parameter posterior uncertainty, and model error.